Taxonomic and Functional Diversity of Benthic Macroinvertebrate Assemblages in Reservoirs of South Korea

Numerous community indices have been developed to quantify the various aspects of communities. However, indices including functional aspects have been less focused on. Here, we examined how community composition varies in response to the environment and discovered the relationship between taxonomic diversity and functional diversity while considering the environment. Macroinvertebrate communities were collected from 20 reservoirs in South Korea. To characterize functional diversity, functional traits in four categories were considered: generation per year, adult lifespan, adult size, and functional feeding groups. Based on their community composition, we classified the reservoirs using hierarchical cluster analysis. Physicochemical and land use variables varied considerably between clusters. Non-metric multidimensional scaling indicated differences between reservoirs and clusters in terms of structure, functional diversity, and environmental variables. A self-organizing map was used to categorize functional traits, and network association analysis was used to unravel relationships between functional traits. Our results support the characteristics of species’ survival strategies such as r- and K-selection. Functional richness exhibited a relationship with taxonomic diversity. Our findings suggest that different types of diversity could play complementary roles in identifying biodiversity. Our findings should prove useful in developing new criteria for assessing freshwater ecosystem health, as well as in evaluating and predicting future alteration of benthic macroinvertebrate communities facing anthropogenic disturbances.


Introduction
To quantify the features of communities and ecosystems, various community indices have been developed. Most indices have traditionally been based on taxonomic richness and abundance, giving rise to taxonomic diversity indices such as the Simpson diversity index [1] and the Shannon diversity index [2]. Taxonomic diversity has the advantage of providing numerical values immediately, without the need for any additional measures other than identification [3]. Based on this advantage, taxonomic diversity has been widely used in various types of biological indices to evaluate ecosystem responses to environmental changes and disturbances [4]. Although taxonomic diversity properly depicts the condition of the systems under consideration, it does not present ecosystem functioning, which is one of two key aspects (i.e., structure and functioning) in ecosystems.
The life history of species reflects their adapted traits in both biological and ecological aspects through natural selection to their environment [5]. Therefore, the functional traits of species display the characteristics of responses of organisms to their environment and are key for understanding community characteristics that are difficult to determine using taxonomic community indices alone [5]. Functional diversity is defined as the diversity of functional traits in a community, and it comprises components of biodiversity that impact how an ecosystem operates or functions [6]. Mason et al. [7] proposed three popular A dataset of benthic macroinvertebrates, which were surveyed at 20 reservoirs in the Nakdong River catchment, South Korea ( Figure 1, Table S1) from 2009 to 2017, was obtained from reports of the Survey on the Environment and Ecosystem of Lakes (SEEL) [25]. The Nakdong River flows into the southern part of the Korean Peninsula and is the longest river in South Korea (total river length: 8567 km), with a catchment area of 30,047 km 2 [26,27]. Macroinvertebrates were collected twice a year (in spring and autumn) using a dredge net (40 cm width, 0.5 mm mesh size) [25]. According to the SEEL protocol [25], samples were collected at two sampling points where water was present throughout the year in each reservoir [28], and the samples were pooled for each reservoir, and average abundance (individuals m 2 /sample/year) data were used in the analyses. The final dataset consisted of 90 benthic macroinvertebrate taxa from 20 reservoirs. abundance (individuals m 2 /sample/year) data were used in the analyses. The final dataset consisted of 90 benthic macroinvertebrate taxa from 20 reservoirs. To evaluate differences in communities based on environmental differences, 19 environmental variables in three categories were used: 12 water quality variables (pH, dissolved oxygen (DO), biochemical oxygen demand (BOD), chemical oxygen demand (COD), total suspended solids (TSS), total nitrogen (TN), total phosphorus (TP), total organic carbon (TOC), water temperature (WT), electric conductivity (EC), chlorophyll-a (Chl-a), and transparency), four geomorphological variables (elevation, storage area, bank height, and circumference) and three land use types (proportions of urban area, agricultural area, and forest). Water quality data were obtained from the Water Environment Information System (WEIS; http://water.nier.go.kr/, accessed on 7 June 2021) of South Korea of the Ministry of Environment in South Korea. Water quality was measured every month at each reservoir according to the WEIS protocol. We used average values of water quality variables in the sampled year in the analyses. Land use categories within a catchment of the study reservoirs were measured from a digital map using QGIS [29]. Physical environmental data were obtained from the Water Management Information System (WAMIS; http://www.wamis.go.kr/, accessed on 16 November 2022).

Functional Traits of Macroinvertebrates
Based on the literature, each taxon of macroinvertebrates was characterized by functional traits in four categories, including the number of generations per year (voltinism), adult lifespan, adult size to maturity, and functional feeding groups of taxa at the genus level (Table 1) [30][31][32][33][34][35]. References used to define functional traits were given in the supplementary information Table S3. The number of generations per year, adult life span, and adult size are related to the population dynamics of the species, whereas functional feeding groups are related to food availability. Voltinism aided in determining the diversity of the age structure, which influences the response patterns of the population to disturbances [36]. Adult life span is related to the reproductive potential of a species [37]. In contrast, adult size incorporates ecological and physiological responses to the local environment and is related to fecundity [38]. To evaluate differences in communities based on environmental differences, 19 environmental variables in three categories were used: 12 water quality variables (pH, dissolved oxygen (DO), biochemical oxygen demand (BOD), chemical oxygen demand (COD), total suspended solids (TSS), total nitrogen (TN), total phosphorus (TP), total organic carbon (TOC), water temperature (WT), electric conductivity (EC), chlorophyll-a (Chl-a), and transparency), four geomorphological variables (elevation, storage area, bank height, and circumference) and three land use types (proportions of urban area, agricultural area, and forest). Water quality data were obtained from the Water Environment Information System (WEIS; http://water.nier.go.kr/, accessed on 7 June 2021) of South Korea of the Ministry of Environment in South Korea. Water quality was measured every month at each reservoir according to the WEIS protocol. We used average values of water quality variables in the sampled year in the analyses. Land use categories within a catchment of the study reservoirs were measured from a digital map using QGIS [29]. Physical environmental data were obtained from the Water Management Information System (WAMIS; http://www.wamis.go.kr/, accessed on 16 November 2022).

Functional Traits of Macroinvertebrates
Based on the literature, each taxon of macroinvertebrates was characterized by functional traits in four categories, including the number of generations per year (voltinism), adult lifespan, adult size to maturity, and functional feeding groups of taxa at the genus level (Table 1) [30][31][32][33][34][35]. References used to define functional traits were given in the supplementary information Table S3. The number of generations per year, adult life span, and adult size are related to the population dynamics of the species, whereas functional feeding groups are related to food availability. Voltinism aided in determining the diversity of the age structure, which influences the response patterns of the population to disturbances [36]. Adult life span is related to the reproductive potential of a species [37]. In contrast, adult size incorporates ecological and physiological responses to the local environment and is related to fecundity [38]. Taxonomic community indices such as taxa richness, abundance, the Shannon diversity index [2], evenness, and the dominance index were calculated for each study reservoir. In addition, functional traits were used to calculate functional diversity such as functional richness (FRic) and functional evenness (FEve) [8,9]. FRic is calculated using the convex hull volume and represents the potential occupancy of niche space [7,8], and FEve is based on the minimum spanning tree and indicates how average species traits are distributed regularly within the occupied trait space [8].

Community Characteristics
To investigate how well the community structure characterizes the reservoirs, a twoway cluster analysis [39] based on community composition was conducted using the Ward linkage method and Bray-Curtis distance. To reduce variation, abundance was log-transformed, and only taxa found at more than two reservoirs (52 taxa) were used in the analyses. One-way analysis of variance (ANOVA) was used to compare differences in environmental conditions among clusters, followed by a post hoc Tukey's HSD multiple comparison test [40]. Non-metric multidimensional scaling (NMDS) was performed using the same data employed in the cluster analysis to determine the relationship between spatial variation in community structure, taxonomic diversity indices, functional diversity indices, and environmental variables [41].

Functional Characteristics of Communities
The communities were analyzed based on the functional traits of species, and the relationships between functional and taxonomic diversity were compared. We used functional traits in four categories (number of generations per year, adult lifespan, adult size, and functional feeding groups), and the analysis included 85 taxa with information of at least three functional trait categories.
To determine the associations between taxa and functional traits, a self-organizing map (SOM) [42] was constructed using 85 taxa data comprising functional traits. The functional traits of each taxon were expressed as a dichotomous transformation (i.e., 0 or 1). SOM output units (40 = 5 × 8) were calculated based on the function 5 × (number of samples) [43,44]. Following the learning process of the SOM, taxa and functional traits were categorized using two-way cluster analysis with the Ward linkage method and the Euclidean distance measure [45].
A network association analysis was performed to evaluate the links between functional attributes. Network association indices such as support, confidence, and lift were used to characterize the association among functional traits. Support for item X was defined as the ratio of transactions in the dataset that contained the item [46]. If items X and Y were present, support (X→Y) refers to the proportions that contained both items X and Y (P (X ∩ Y)). Confidence (X→Y) denotes the conditional probability that Y will be adopted if X is adopted (P (X ∩ Y)/P(X)) [46]. The maximum confidence level was 1, indicating that Y was always adopted when X was adopted. Lift (X→Y) is the ratio of the probability of choosing Y when choosing X to the probability of choosing Y (P (X ∩ Y)/P(X) P(Y)) [47]. This indicates the strength of the link between X and Y. Rules with one-on-one correspondence between functional traits were used with 0.1 of minimum support and 0.8 of minimum confidence as default as well as 0.05 of minimum support and 0.5 of minimum confidence to include more associated networks.
The Spearman rank correlation coefficient was calculated to evaluate the relationship between taxonomic and functional diversity indices. Furthermore, a generalized additive model (GAM) was employed to characterize changes in taxonomic and functional diversity as a function of community indices by accounting for the gradient of environmental variables [48].
All statistical analyses were conducted using the R software (https://www.r-project. org/, accessed on 23 June 2022) with the following related packages: functional diversity indices with the dbFD function in the FD package [49,50], hierarchical cluster analysis and NMDS with the vegan package [51], SOM with the kohonen package [52], ANOVA with the stat package [53], Tukey's HSD multiple comparison test with the PMCMRplus package [54], network association analysis with the arules package [55], and GAM with the mgcv package [56,57].
Based on the similarity of their community compositions, the 20 reservoirs were divided into three clusters (S1-S3; Figure 2). Furthermore, the 52 taxa identified in the dataset were classified into four groups (T1-T4) based on similarities in their abundance and occurrence in the study reservoirs. The dominant taxon in all clusters was Chironomus spp.; however, the second dominant taxa varied by cluster, such as Limnodrilus spp., Micronecta spp., and P. paucidens in clusters S1 through S3, respectively. In group T4, cluster S1 (BnSn and JnCh) had greater taxonomic richness and abundance. The taxa in group T2 were widely distributed across the study reservoirs. Taxonomic abundance was low in T1 and T3.
Physicochemical environmental conditions varied across the reservoirs ( Figure 3). Among water quality variables, BOD, COD, TP and TOC were the highest in cluster S1 (ANOVA, p < 0.05), whereas TN was the highest in cluster S2 (ANOVA, p < 0.05). Regarding land use, the proportion of forest was the lowest and that of the agricultural area was the highest in cluster S1 (ANOVA, p < 0.05). Cluster S1 had higher EC and lower transparency, elevations and bank heights than the other clusters, although these differences were not significant (ANOVA, p = 0.14 for both elevations and bank heights).
third most dominant taxa in each cluster characterized the differences in the clusters (Figure 4b). Axis 1 of the NMDS was positively correlated with dominance and FEve, and negatively correlated with richness, abundance, Shannon diversity, and FRic ( Figure 4c). However, FEve and dominance had positive and negative correlations along NMDS axis 2, respectively. Nutrient conditions (BOD, COD, TP and TOC) and the proportion of agricultural area were higher in cluster S1 than in the other clusters, whereas the proportion of forest and elevation were lower in cluster S1 ( Figure 5). Moreover, elevation and the proportion of forest were high in clusters S2 and S3.   The NMDS ordination reflected differences in the communities and the environmental conditions of the reservoirs (Figures 4 and 5). Cluster S1, which had a distinct community composition, was isolated from the other clusters (Figure 4a). The second and the third most dominant taxa in each cluster characterized the differences in the clusters (Figure 4b). Axis 1 of the NMDS was positively correlated with dominance and FEve, and negatively correlated with richness, abundance, Shannon diversity, and FRic ( Figure 4c). However, FEve and dominance had positive and negative correlations along NMDS axis 2, respectively. Nutrient conditions (BOD, COD, TP and TOC) and the proportion of agricultural area were higher in cluster S1 than in the other clusters, whereas the proportion of forest and elevation were lower in cluster S1 ( Figure 5). Moreover, elevation and the proportion of forest were high in clusters S2 and S3.

Functional Traits and Diversity
The 85 taxa were classified into four groups (t1-t4) through the SOM learning process based on similarities in their functional trait composition ( Figure 6). More than half of Hemiptera taxa were in group t1. Non-Insecta taxa, Diptera, Coleoptera and Odonata were typically included in groups t1 and t4, whereas Ephemeroptera taxa were in groups t2 and t3.

Functional Traits and Diversity
The 85 taxa were classified into four groups (t1-t4) through the SOM learning process based on similarities in their functional trait composition ( Figure 6). More than half of Hemiptera taxa were in group t1. Non-Insecta taxa, Diptera, Coleoptera and Odonata were typically included in groups t1 and t4, whereas Ephemeroptera taxa were in groups t2 and t3.  Table 1.
Based on the weight vectors of the trained SOM, functional traits were classified into three clusters (F1-F3; Figure 6b). In F1, taxa in groups t2 to t4 exhibited a highly frequent occurrence of functional traits similar to those in V2; in F2, infrequent functional traits were included similar to those in F1, F5 and S2; in F3, taxa in group t1 and t4 had a high frequency of functional traits such as S3.
Group t1 was characterized by large body sizes (S3; Figure 6b), and the functional feeding groups differed between taxa. Group t2 consisted of taxa with small body sizes (S1) and short lifespans (A1), as well as one generation per year (V2). Group t3 was composed of gatherer-collectors (F2) with one or more generations (V2-V3) per year and medium lifespans (A2). Patterns in group t4 with one generation per year (V2) were distinguished by large body sizes (S3), as well as various feeding habits.
Network association analysis revealed relationships between functional traits (Table 2, Figure 7). Three distinct associations were identified with minimum support of 0.1 and minimum confidence of 0.8. Herbivorous feeding pattern (F3) resulted in a large body size (S3). Less than one generation per year (V1) also associated with large body size (S3) and related to long lifespans (A3; Figure 7a). The network association results supported those of the hierarchical cluster analysis (Figure 6b), especially functional traits in group F3. Moreover, four modules with modularity of 0.32 were generated with minimum support of 0.05 and minimum confidence of 0.5 (Figure 7b). Module M1 showed patterns similar to those of the functional trait cluster analysis in group F3 (Figure 6b). Taxa with fewer  Table 1.
Based on the weight vectors of the trained SOM, functional traits were classified into three clusters (F1-F3; Figure 6b). In F1, taxa in groups t2 to t4 exhibited a highly frequent occurrence of functional traits similar to those in V2; in F2, infrequent functional traits were included similar to those in F1, F5 and S2; in F3, taxa in group t1 and t4 had a high frequency of functional traits such as S3.
Group t1 was characterized by large body sizes (S3; Figure 6b), and the functional feeding groups differed between taxa. Group t2 consisted of taxa with small body sizes (S1) and short lifespans (A1), as well as one generation per year (V2). Group t3 was composed of gatherer-collectors (F2) with one or more generations (V2-V3) per year and medium lifespans (A2). Patterns in group t4 with one generation per year (V2) were distinguished by large body sizes (S3), as well as various feeding habits.
Network association analysis revealed relationships between functional traits (Table 2, Figure 7). Three distinct associations were identified with minimum support of 0.1 and minimum confidence of 0.8. Herbivorous feeding pattern (F3) resulted in a large body size (S3). Less than one generation per year (V1) also associated with large body size (S3) and related to long lifespans (A3; Figure 7a). The network association results supported those of the hierarchical cluster analysis (Figure 6b), especially functional traits in group F3. Moreover, four modules with modularity of 0.32 were generated with minimum support of 0.05 and minimum confidence of 0.5 (Figure 7b). Module M1 showed patterns similar to those of the functional trait cluster analysis in group F3 (Figure 6b). Taxa with fewer than one generation per year (V1) were associated with predatory feeding patterns (F4) and long lifespans (A3). Herbivores (F3) had a large body size (S3), whereas small body size (S1) resulted in gatherer-collectors (F2) and short lifespans (A1). The characteristics of more than one generation per year (V3) and gatherer-collector feeding habits (F2) were linked. Moderate characteristics, shredders (F5), and filterer-gatherers (F1) were included in modules M3 and M4, respectively. Members of module M4 were moderate characteristics. Table 2. Association rules for functional traits in Figure 7a. Abbreviations of traits were described in Table 1.

No.
Antecedent Consequent Lift Support (%) than one generation per year (V1) were associated with predatory feeding patterns (F4) and long lifespans (A3). Herbivores (F3) had a large body size (S3), whereas small body size (S1) resulted in gatherer-collectors (F2) and short lifespans (A1). The characteristics of more than one generation per year (V3) and gatherer-collector feeding habits (F2) were linked. Moderate characteristics, shredders (F5), and filterer-gatherers (F1) were included in modules M3 and M4, respectively. Members of module M4 were moderate characteristics.  Table 1. Table 2. Association rules for functional traits in Figure 7a. Abbreviations of traits were described in Table 1
With the indices demonstrating significant correlations (Figure 8), the relationships between taxonomic and functional diversity indices were analyzed further by taking environmental factors into account (Figure 9). Richness and abundance explained 52.0% and 60.3% of the deviance in FRic by GAM, respectively. Environmental variables tended to increase or decrease as taxa richness, abundance and FRic increased. Taxa richness, abundance and FRic were positively correlated with BOD. TP was positively related to taxa richness and FRic, whereas these indices negatively correlated with proportion of forest and elevation. With the indices demonstrating significant correlations (Figure 8), the relationships between taxonomic and functional diversity indices were analyzed further by taking environmental factors into account (Figure 9). Richness and abundance explained 52.0% and 60.3% of the deviance in FRic by GAM, respectively. Environmental variables tended to increase or decrease as taxa richness, abundance and FRic increased. Taxa richness, abundance and FRic were positively correlated with BOD. TP was positively related to taxa richness and FRic, whereas these indices negatively correlated with proportion of forest and elevation.   With the indices demonstrating significant correlations (Figure 8), the relationships between taxonomic and functional diversity indices were analyzed further by taking environmental factors into account (Figure 9). Richness and abundance explained 52.0% and 60.3% of the deviance in FRic by GAM, respectively. Environmental variables tended to increase or decrease as taxa richness, abundance and FRic increased. Taxa richness, abundance and FRic were positively correlated with BOD. TP was positively related to taxa richness and FRic, whereas these indices negatively correlated with proportion of forest and elevation. Figure 9. Relationship between taxonomic (taxa richness (a), abundance (b)) and functional diversity and environmental variables. Only statistically significant relationships are plotted. GAM was Figure 9. Relationship between taxonomic (taxa richness (a), abundance (b)) and functional diversity and environmental variables. Only statistically significant relationships are plotted. GAM was applied to analyze the relationship between taxonomic and functional diversity (black lines) with 95% confidence intervals (grey-shaded areas). Point color presents relative value of each environmental variable (light: low value and dark: high value).

Discussion
We evaluated the association between the taxonomic and functional diversity of benthic macroinvertebrates with environmental variables in 20 reservoirs. By linking environmental conditions, we were able to identify changes in community structure in response to environmental changes. Linkages between functional features and the association between biodiversity indicators were also identified.

Taxonomic Community Structure
Water quality has a significant influence on the composition of taxonomic communities. In particular, groups were separated by the parameters BOD, COD, TN, TP and TOC (Figure 3). Based on the water quality standard of South Korea (https:/water.nier.go.kr/, accessed on 7 June 2021), differences in COD, TN, TP and TOC standards existed among the study reservoirs. In this study, the BnSn and JnCh reservoirs (S1) had a distinct community composition as compared to the other reservoirs ( Figure 2). Cluster S1 (BnSn and JnCh) also had a relatively higher BOD, with an average of 3 mg/L as compared to the other clusters (whose BOD was approximately 2 mg/L), highlighting the effects of BOD on community composition (Figure 3).
There were numerical disparities in EC and transparency across clusters (Figure 3), which influenced community structure in various ways. EC indicates the overall water quality by integrating the total number of dissolved ions in water. On the other hand, the transparency of water influences food web dynamics in aquatic ecosystems [58] and functional traits communities, such as foraging strategies and success [59].
TN, TP and TOC are known to be closely related to community composition [60]. Nitrogen and phosphorus enable tolerant species to thrive in lakes [61]; in particular, TP regulates the taxonomical and functional beta diversity of macroinvertebrates [60]. In addition, TOC is influenced by inputs from the watershed [62]. Water quality is closely related to land use and cover, in which the proportions of the urban landscape, agricultural land, and forest areas are particularly important [63,64]. In our study, we found that the landscape had a distinct composition in S1, where the community differed from those in the other clusters (Figure 3). In addition, landscape composition showed similar or opposite tendencies to water quality variables such as BOD, COD, and TOC ( Figure 5).
The ecosystem health index based on benthic macroinvertebrates is closely related to chemical and physical gradients and displays a strong correlation with land use, suggesting that land use indirectly induces considerable changes in macroinvertebrate communities [65]. In inland wetlands, the main drivers of water quality degradation are nutrient and sediment inputs, which are associated with altered land use [66]. Human-induced land use has been the major source of nitrogen and phosphorus in lakes, while forests served to reduce nitrogen and phosphorus contents in lakes [67].

Relationships between Functional Traits of Macroinvertebrates
Functional traits were classified by the SOM learning process based on the presence or absence of traits in the observed species, and their association was analyzed using network association analysis (Figures 6b and 7). The taxonomic system reflects biological characteristics in the community [32], whereas the functional traits of species, such as the combination of size, reproductive traits (r to K strategies), and food and feeding habitats, could explain taxonomic groups because functional traits reflect the evolutionary responses of species to their environmental changes [68]. Some functional traits have similarities or trade-offs, such as body size and fecundity [38].
Even though the taxonomic system and biological characteristics are weakly related [32,68], the classification of functional traits was not completely consistent with that of the taxonomic characteristics ( Figure 6), indicating that functional and taxonomic community functions differ. Taxa that survive environmental filtering have several linked traits that allow them to adapt well to their environment [69]. Our findings were consistent with the characteristics of rand K-selection taxa [70]. Rapid development and early reproduction, and short lifespan are related to r-selection, whereas K-selection taxa are related to slower development, delayed reproduction, large body size, and long lifespan. High voltinism is caused by a combination of fast development and rapid reproduction.
There is a trade-off between body size and voltinism [71], which is reflected in the findings of our study (Figures 6b and 7). Although Zeuss et al. [71] found that climate could influence this relationship, our study sites had a comparable climate; therefore, a minor temperature variation may not have affected this relationship in this study. The close size-voltinism relationship is usually valid in areas with similar climates. Species with frequent voltinism grow quickly, and the energy necessary for rapid growth lowers the priority of other abilities [72].
Body size and trophic level are closely linked in aquatic ecosystems [73]. In this study, large body size was closely related to predatory feeding habits, whereas small body size was related to gatherer-collector feeding habits (Figure 7). Gatherer-collectors, which feed on organic matter or primary producers, have become the primary food source for predatory insects [74], and have a low trophic level. Large species require more energy than small species since they require more energy to grow. Predators tend to have a larger body size and require a high-caloric energy source (food). As the size of the predator increases, prey size tends to increase [75]; the consumption rate also increases with the predator body mass [76]. Therefore, predator size appears to be directly related to food consumption.
Taxa with high trophic levels and predatory feeding habits are density-dependent on prey in that food energy efficiency is high, but quantity is limited [77]. Density dependence is observed not only in interspecific competition but also in intraspecific competition. There are limitations, however, because size-dependent intraspecific predation occurs between individual predators with low voltinism; therefore, larvae have varied size structures [78]. Furthermore, the larger the body size, the smaller the intraspecific adult size variation, implying that a large-sized species requires a high-quality habitat in which body size may grow; hence, the suitability of the environment could be evaluated by the presence of large predators [79].
Movement capability is related to tolerance to environmental changes [80], and tolerance is associated with life span. A taxon with a reduced lifespan is vulnerable to environmental changes [80], seemingly due to difficulties in surviving for a sufficient period of time until a suitable environment is built; therefore, they are replaced by species with greater fitness. In consideration of the characteristics of r-selection species with short life spans, r-selection species with a higher number of voltinisms are hypertrophied in a fit environment, or do not fit and occur less frequently. Therefore, in the absence of major disturbances, a community with longer-lived species shows less variation in community structure and decreases the temporal heterogeneity of the community [81]. Considering the slow growth and long lifespans of K-selection species, their presence at moderate perturbation aids in lowering community heterogeneity, since slow-growing species are more readily adapted to disturbances [82]. As such, functional traits also influence temporal changes in community structure.

Taxonomic Diversity and Functional Diversity
Among the correlations between taxonomic and functional diversity indices, FRic was related to taxonomic indices, whereas FEve was not (Figure 8). Taxonomic and functional diversity have a positive relationship, although the two are partially independent [83]. Although Hacala et al. [84] showed that taxonomic and functional diversity were strongly positively correlated, they are based on fundamentally different perspectives, despite of some overlap. Devictor et al. [16] mentioned the mismatching role of taxonomic and functional diversity to ecosystem estimation; that is, functional underestimation and taxonomic overestimation in ecoregions. Therefore, while taxonomic diversity indices cannot be completely replaced by functional diversity indices [85], they can be used to examine ecosystems.
Although there were few statistically significant linear relationships between the functional diversity indices and environmental variables in this study (Figure 9), there has been a lot of research on this topic in recent years, and some studies have found that functional diversity is partially affected by the environment. For example, functional richness declines with disturbances such as urbanization [86]. The functional diversity index fluctuates in response to changes in the habitat environment in addition to disturbance, having a higher impact than spatial vectors [87]. FRic is affected by water temperature, substratum, precipitation, and water quality parameters such as pH and total nitrogen [87,88]; FEve is also affected by substratum [88]. Furthermore, partial functional diversity such as FEve is occasionally independent of disturbances such as land-use changes [89]. Functional diversity did not exhibit a linear relationship with changes in the environment in this study; however, it was expected to feature other forms of relation, such as unimodal.

Limitations of the Study
We revealed the relationships between taxonomic and functional diversities and environmental variables in reservoirs. However, we are aware of the limitations of this study which are commonly observed in studies using datasets obtained from a public database or literature. Our dataset was from the SEEL. In this monitoring program, benthic macroinvertebrates were collected with a dredge net by drawing the sediment on the bottom of reservoirs and to prevent loss of active macroinvertebrates [28,90]. Kick net sampling may be more efficient at sampling benthic macroinvertebrates than other methods, but kick net sampling is only feasible in shallow water [90]. Sampling with a dredge is typically used to characterize benthic macroinvertebrate assemblages in lakes [91]. de Mendoza and Catalan [92] determined the sampling points of lakes according to the habitat type, and integrated samples in each lake for further analyses, supporting the sampling strategy of our study. Thus, our sampling of benthic macroinvertebrates in the reservoirs with a dredge is justifiable [28,93]. Despite these potential limitations, our results are comparable between taxonomic and functional diversities and environmental variables in the study reservoirs.

Conclusions
We discovered the relationships between taxonomic and functional diversity and environmental variables in reservoirs. Communities were classified into three clusters based on taxonomic diversity, reflecting differences in environmental conditions as well as taxonomic and functional diversity. Functional traits in the four categories were classified into three groups by the SOM learning process, demonstrating a connection, similarities, or a trade-off relationship between functional traits such as the number of generations per year, adult size, and lifespan. Network association analysis revealed relationships between functional traits, such as a medium lifespan with one generation per year. Our findings support the characteristics of species' survival strategies such as r-and K-selection. Among the functional diversity indices, FRic had the strongest relationship with taxonomic diversity. When environmental variables were considered, a linear relationship between both taxonomic diversity and functional diversity indices such as taxa richness, abundance and FRic and environmental conditions was discovered. Taxonomic and functional diversity showed both comparable and different trends, indicating that both could be used to complement ecosystem analysis. From this perspective, our findings should prove useful in developing new criteria for assessing freshwater ecosystem health, as well as in evaluating and predicting future alteration of benthic macroinvertebrate communities facing anthropogenic disturbances.