Effects of Land Use Types on Community Structure Patterns of Benthic Macroinvertebrates in Streams of Urban Areas in the South of the Korea Peninsula

Benthic macroinvertebrates were collected from streams located in an urban area from regions featuring different environmental conditions. Physicochemical variables and land use types pertaining to sampling sites were analyzed concurrently. Multivariate analyses (cluster analysis and non-metric multidimensional scaling) and rank-abundance diagrams were used to characterize community patterns to assess ecological integrity in response to environmental conditions. Species composition patterns were mainly influenced by both the gradient of physicochemical variables (e.g., altitude, slope, conductivity) and the proportion of forest area. Community structure patterns were further correlated to the proportion of urbanization and to biological indices (e.g., diversity, number of species). Land use preferences of benthic species were identified based on the indicator values and weighted averaging regression models. Plecoptera species were representative of undisturbed streams in forest areas, whereas Tubificidae species and filtering collector caddis flies were indicator taxa in severely polluted and agricultural areas, respectively. The analyses of community structures and indicator species effectively characterized community properties and ecological integrity following natural and anthropogenic variability in urban stream ecosystems.


Introduction
Due to unprecedented industrial development and human aggregation, ecosystems have been severely damaged and destabilized both locally and globally.Urbanization causes rapid changes in land use via expansion of commercial/residential/industrial areas along with decimation of forests and natural areas [1].The anthropogenic impact of urbanization upon aquatic ecosystems is critical, considering the importance of resource usage and the diffusion of pollutants to the surrounding environment [2,3].Ecological degradation of stream drainage was especially affected by urbanization, causing the so-called 'urban stream syndrome' [3].The symptoms of the urban stream syndrome have been broadly characterized in various fields, including physicochemical and ecological processes such as hydrology [3,4], biodiversity conservation [5][6][7], and ecosystem processes [8,9].
The conversion of natural areas (e.g., forest) to urban or agricultural areas severely affects both in-stream habitat and macroinvertebrate communities in various ways [10,11].Both the abundance Water 2016, 8, 187 2 of 18 and the distribution of aquatic organisms are strongly influenced by land use [12,13].Benthic macroinvertebrates adapt to the environmental changes caused by various anthropogenic impacts, resulting in an increase of tolerant species in the community.Therefore, community composition serves as a good predictor of habitat quality and biotic integrity [14,15].During the last decade, studies demonstrated that a change in land use was one of the main driving forces behind the loss of regional biodiversity [11,16].
Community data are difficult to analyze because communities consist of numerous species varying in a complex and stochastic manner in response to environmental factors.Computational methods have been developed to analyze community data efficiently; these include predictive models [17], multivariate statistical methods [18], statistical learning models [19,20], and species abundance models [21].Many studies have performed multivariate analyses on freshwater benthic community data to assess water quality [18,22], spatial and temporal dynamics [23,24], trait patterns [25], and community response to disturbance [18,26].A rank abundance diagram (RAD), a type of species abundance model, has been applied to community data to elucidate alterations in community structure resulting from environmental change [27].Several studies using RADs have been applied to various taxa, including freshwater macroinvertebrates [26,28].
In this study, we extended the results of previous studies [21,28] to further investigate the changes in macroinvertebrate community patterns in response to urbanization in the southeast part of the Korean Peninsula.Specifically, we presented the environmental and biological characteristics, and examined the effects of anthropogenic disturbance on benthic macroinvertebrate communities by concurrently measuring hydrological/physical/chemical factors and land use types in various streams located in the Busan metropolitan area.

Field Sampling
Benthic macroinvertebrates were collected at 46 sampling sites in the Busan metropolitan area, located in the southeast part of the Korean Peninsula (Figure 1a).The geography of the Busan metropolitan area is covered in mountainous regions (47%), and the southeast part of Korea features a temperate climate with high rainfall in July because of the East Asian Monsoon.Sampling was carried out twice, namely, in late summer (August) and in winter (November-December), in 2013 and 2014.Sampling sites were chosen across various altitudes (1-315 m) and for different land use types.Among the 46 sampling sites, 10 were located on four tributaries of the Nakdong River, which flows through the Busan metropolitan area, and is the longest river in South Korea.Twenty-one sampling sites were located in the Suyeong River of mid-eastern Busan, and 15 sites were located on short streams flowing directly into the East Sea (Figure 1).The overall environmental conditions of each sampling site are shown in Appendix.
We collected macroinvertebrate samples in triplicate from a riffle habitat within 10 m reaches using a Surber sampler (30 ˆ30 cm 2 , 500 µm mesh).Ten environmental variables, including geological, landscape, hydrological, and physicochemical factors, were measured concurrently (Appendix).Water quality parameters, including electrical conductivity (YSI 30, Yellow Springs Instruments, Yellow Springs, OH, USA) and dissolved oxygen (DO) (YSI 550A, Yellow Springs Instruments, Yellow Springs, OH, USA), were measured in situ, along with hydrological variables such as velocity (Swoffer 2100 LX, Swoffer Instruments, Seattle, WA, USA), depth, and substrata.Hydrological variables were measured for each replicate.Substrates collected by the Surber sampler were classified in five categories according to Jowett et al. [29]; a substrate index was calculated based on the substrate composition [30].Topological and geological variables, including altitude, slope, and land use types, were extracted from a digital map obtained from the Korea Ministry of Environment using ArcGIS version 9.3 (ESRI).Proportions (%) of land use types (specifically, forest, urban, agriculture, grassland, wetland, and bare land) within a 1 km-long riparian zone (100 m width) were calculated from the digital map.calculated from the digital map.This riparian sub-corridor was selected, because variation in land use is often greater at catchment scales than it is at riparian scales [2,31].
Benthic macroinvertebrate specimens were mostly identified either to the species level or to the lowest possible taxonomic level according to Han et al. [32], Won et al. [33], and Yoon [34].In addition, Merritt and Cummins [35], and Brigham et al. [36] were used as ancillary literature.Chironomidae and Oligochaeta were identified to the family or class level, owing to difficulty in identification.

Data Analysis
The community data pooled at each site were used in the analyses.The Shannon diversity index, using a common logarithm [37], the dominance index [38], and Ephemeroptera, Plecoptera, and Trichoptera (EPT) richness, were calculated at each sampling site from the community data.Spearman rank correlation coefficients were obtained to show associations between the biological indices and the environmental variables.
Considering the complexity of the community data due to environmental variability, we applied multivariate analysis and RAD to determine community patterns.Before applying cluster analysis and nonmetric multidimensional scaling (NMS), the community abundance data were log transformed with ln(x + 1) to reduce the high variance caused by the exceptionally large number of individuals collected at the sampling sites, and to fulfill the assumption of normality for statistical analyses.The sampling sites were classified with a hierarchical cluster analysis, using Ward's linkage method with a Bray-Curtis distance measure, based on the similarity of their species compositions.A Multi-response Permutation Procedure (MRPP) was used to test whether or not there were statistically significant differences among the clusters [39].NMS was then applied to characterize the overall relationships among communities in association with specific environmental factors.Cluster analyses and NMS were performed using the PC-ORD software [40].A Monte Carlo test was implemented using both randomized data (one thousand runs) and real data to estimate the significance of the output of the NMS [40].Kruskal-Wallis test and a multiple comparison test (Dunn's post hoc test) were applied to identify differences among clusters using the environmental variables and biological water quality indices (p < 0.05).
Slope lines of RADs were obtained to show community responsiveness to environmental conditions [26,41,42].The relative abundance of species (i) in relation to the total abundance for all species in each sampling site was obtained and rearranged according to rank (in order from highest to lowest proportion) to calculate RAD slopes.The slope (S) was calculated with Equation (1): where n is the number of species and pi is the relative abundance of species i according to the rank.Diversity is maximized when the slope of the rank abundance curve approaches zero (species having a greater chance of being evenly represented), while slopes become increasingly negative when total Benthic macroinvertebrate specimens were mostly identified either to the species level or to the lowest possible taxonomic level according to Han et al. [32], Won et al. [33], and Yoon [34].In addition, Merritt and Cummins [35], and Brigham et al. [36] were used as ancillary literature.Chironomidae and Oligochaeta were identified to the family or class level, owing to difficulty in identification.A list of taxa identified in the study sites are given with their literature in the Supplementary Material (Table S1).

Data Analysis
The community data pooled at each site were used in the analyses.The Shannon diversity index, using a common logarithm [37], the dominance index [38], and Ephemeroptera, Plecoptera, and Trichoptera (EPT) richness, were calculated at each sampling site from the community data.Spearman rank correlation coefficients were obtained to show associations between the biological indices and the environmental variables.
Considering the complexity of the community data due to environmental variability, we applied multivariate analysis and RAD to determine community patterns.Before applying cluster analysis and nonmetric multidimensional scaling (NMS), the community abundance data were log transformed with ln(x + 1) to reduce the high variance caused by the exceptionally large number of individuals collected at the sampling sites, and to fulfill the assumption of normality for statistical analyses.The sampling sites were classified with a hierarchical cluster analysis, using Ward's linkage method with a Bray-Curtis distance measure, based on the similarity of their species compositions.A Multi-response Permutation Procedure (MRPP) was used to test whether or not there were statistically significant differences among the clusters [39].NMS was then applied to characterize the overall relationships among communities in association with specific environmental factors.Cluster analyses and NMS were performed using the PC-ORD software [40].A Monte Carlo test was implemented using both randomized data (one thousand runs) and real data to estimate the significance of the output of the NMS [40].Kruskal-Wallis test and a multiple comparison test (Dunn's post hoc test) were applied to identify differences among clusters using the environmental variables and biological water quality indices (p < 0.05).
Slope lines of RADs were obtained to show community responsiveness to environmental conditions [26,41,42].The relative abundance of species (i) in relation to the total abundance for all species in each sampling site was obtained and rearranged according to rank (in order from highest to lowest proportion) to calculate RAD slopes.The slope (S) was calculated with Equation (1): where n is the number of species and p i is the relative abundance of species i according to the rank.Diversity is maximized when the slope of the rank abundance curve approaches zero (species having a greater chance of being evenly represented), while slopes become increasingly negative when total diversity decreases [26].For convenience of expression, we used the absolute value to present the slopes.Therefore, larger values represented steeper rank abundance curves, thus showing a tendency toward less diverse assemblages.
The indicator species analysis was conducted to define representative species in different land use types.Indicator values (IndVal) for each species i in cluster k were calculated based on Equation (2) [43,44]: where Nindividuals ki is the mean abundance of species i across the sites pertaining to cluster k and Nindividual s+k is the sum of the mean abundances of species i within the various clusters.Nsites kj is the number of sites in cluster k where species i is present and Nsites k+ is the total number of sites in that cluster.IndVal ranges from 0 (no indication) to 100 (perfect indication).A Monte Carlo method using 300 permutation runs was also implemented to evaluate the statistical significance of IndVal.
Weighted averaging regression models (WARMs) were applied to quantify the preference of each indicator species for each land use type [44,45].The optimal preference value of the species was calculated as the mean of the measured land use proportion weighted by the abundance of this species at the study sites, according to Equation (3) [46]: where WA i is the weighted average (estimate of optimum preference of species i) of each land use type, y j is the proportion of each land use type at site j, and x ij is the abundance of species i at site j.The tolerant condition (TOL i ) of each land use type for species i was calculated as the weighted standard deviation of the species abundance at all sites according to Equation (4) [46]: WARMs were performed using the C2 program environment [47].The coefficients of determination (R 2 ) were used to estimate the precision of WARMs.Model errors were estimated by bootstrapping with 1000 cycles [48].If the regression coefficients estimated by the bootstrapping model, and by the original model, were similar (about 80%), then the results of WARMs could be used [45].

Community Composition
A total of 136 taxa were identified to the species or genus level, except for the family Chironomidae and class Oligochaeta, with an average of 22.2 (˘10.5)taxa, and an average density of 1386.7 individuals/m 2 (˘1380.7)at each site.Twenty-six non-insecta species (19%) were collected, including Physa acuta, Asellus sp., and Semisulcospira libertina.Three species (Baetis fuscatus, Hydropsyche kozhantschikovi, and Physa acuta), and the species of family Chironomidae and class Oligochaeta were widely distributed, and showed a high frequency of occurrence (>70%) in the study sites.The Ephemeroptera species, Baetis fuscatus and Baetiella tuberculata, as well as the Trichoptera species Hydropsyche kozhantschikovi were the most abundant.More than half of the species had a low frequency of occurrence (<10%), and ten singletons (7%) were collected, including four Trichoptera species.The values of community indices, including the number of species and diversity, were higher for the sampling sites at mountain streams than in the other areas.Biological indices varied according to disturbance effects (e.g., high values of conductivity and proportion of urban land use) of sampling sites.

Relationships between Variables
Altitude and slope had strong negative correlations with electrical conductivity (r = ´0.69 and ´0.52 respectively; p < 0.01), but positive correlations with SI (r = 0.61, p < 0.01 and r = 0.31, p < 0.05, respectively) (Table 1).Low altitude and slope indicated the strong anthropogenic impacts.Conductivity, representing the anthropogenic impact, was negatively correlated with SI (r = ´0.57;p < 0.01) and velocity (r = ´0.50;p < 0.01).The proportion of forest area in the riparian zone was positively correlated with the variables related to the mountain streams (altitude: r = 0.70, slope: r = 0.71 with p < 0.01).The overall proportion of urban area had a negative correlation with natural environmental factors, including altitude (r = ´0.62;p < 0.01) and slope (r = ´0.59;p < 0.05), but it was positively correlated with anthropogenic factors, such as conductivity (r = 0.58; p < 0.01).
Conductivity, however, was consistently and strongly negatively correlated with most biological indices (p < 0.01), including the number of species (r = ´0.65),diversity (r = ´0.56),dominance (r = ´0.52),and EPT% (r = ´0.78).Depth and conductivity were similarly correlated with the biological indices.SI and velocity, however, were negatively correlated with depth.The proportion of forest area was strongly correlated with EPT% (r = 0.52; p < 0.01).Negative correlations between urban and biological indices were detected in all cases except for the dominance index.No strong correlation was observed between agricultural land use and the biological variables (Table 1).

Community Classification
The communities were classified into four clusters based on Bray-Curtis dissimilarities in the dendrograms of hierarchical cluster analysis (MRPP, A = 0.42, p < 0.01) (Figures 1 and 2).The classification reflected community similarities as well as environmental gradients.Altitude, slope, and DO were considerably higher in cluster 1, indicating that communities in this cluster were mainly collected from mountain streams in forest-dominant areas (Table 2).In contrast, communities in cluster 4 were from lowland areas with urban environments, having low values for altitude, slope, and DO.The values of the environmental variables in clusters 2 and 3 fell in the middle range between clusters 1 and 4. Conductivity was significantly lower in cluster 1 than in cluster 4. SI and the current velocity were also significantly higher in cluster 1 than in the other clusters, while depth was notably high in cluster 4. Clusters 2 and 3 predominated in the agricultural areas in riparian zones (Table 2).It is noteworthy that diversity values in the highest range were observed in cluster 2 rather than cluster 1.Overall, the samples in cluster 1 from the mountainous areas had the highest values for the biological parameters including density, number of species and EPT%, whereas, except for the dominance index, the lowest values for the biological parameters were observed in cluster 4 (Table 2).

Community Ordination
For the NMS ordination, the sampling sites were distinguished based on similarities in species composition of communities (Figure 3).The final stress value obtained in the NMS ordination was 18.9 (two axes), which is an acceptable stress value for ecological community data, since it falls within the range of 10-20.A two-dimensional ordination explained 66.7% of the variance (R 2 = 0.48 for axis 1, R 2 = 0.19 for axis 2), and the Monte Carlo test implied that the axes explained significantly more variance than by chance alone (p = 0.01).The sampling sites in cluster 1 were aggregated on the upper left part of map, whereas sampling sites in cluster 4 were located on the bottom right part.Clusters 2 and 3 lay on the middle part, between clusters 1 and 4, with cluster 2 to the left and cluster 3 to the right.
The environmental factors had significant correlations with the NMS axes (Figure 3).The main variability in the community data was observed along axis 1 in response to anthropogenic (urban area) and habitat factors.The environmental factors including slope, conductivity, altitude, and the proportion of forest area were related to the bottom and top positions along axis 2. The variability of axis 1 was positively correlated with the proportion of urban area (r = 0.50; p < 0.05), depth (r = 0.44; p < 0.05), and conductivity (r = 0.26; p < 0.05) on the right side with samples in cluster 4, but was negatively correlated with the substrate index (r = ´0.52;p < 0.05) and the current velocity (r = ´0.42;p < 0.05) on the left side of NMS ordination.The variability on axis 2 was lesser than that on axis 1, but it was positively correlated with the slope (r = 0.69; p < 0.05), followed by the proportion of forest area (r = 0.61; p < 0.05), altitude (r = 0.43, p < 0.05), and substrate size (r = 0.41; p < 0.05) on the upper side with samples in cluster 1, and negatively correlated with conductivity (r = ´0.72;p < 0.05) and the proportion of the urban area (r = ´0.44;p < 0.05) on the lower side (Figure 3).The proportion of

Species Abundance Distribution
RADs were obtained as log abundance scales across different clusters (Figure 4), and reflected community conditions pertaining to specific clusters.RADs have low slope and long length, indicating even, highly diverse communities.All slopes were fitted with linear models using higher values of R 2 (average of 0.90 ± 0.07).The slopes were similarly lower in the sampling sites of clusters 1 and 2, which represented the less polluted streams.RADs at the polluted sites in cluster 3 had relatively steeper slopes, with fewer species compared with the less polluted sites in clusters 1 and 2 (Figures 4 and 5).Communities in cluster 4, representing severely polluted sites, had steeper RAD curves and a minimal number of species.All slope values ranged from −0.12 (at YBV) to −1.90 (at YPY).The absolute values of the slopes in clusters 1 and 2, representing the less polluted sites, fell within the range of 0.12-0.32,whereas the values in cluster 3 were relatively higher (within the range of 0.20-1.00).The slopes in cluster 4 were significantly higher than those of clusters 1 and 2, ranging from 0.64-1.90.

Species Abundance Distribution
RADs were obtained as log abundance scales across different clusters (Figure 4), and reflected community conditions pertaining to specific clusters.RADs have low slope and long length, indicating even, highly diverse communities.All slopes were fitted with linear models using higher values of R 2 (average of 0.90 ˘0.07).The slopes were similarly lower in the sampling sites of clusters 1 and 2, which represented the less polluted streams.RADs at the polluted sites in cluster 3 had relatively steeper slopes, with fewer species compared with the less polluted sites in clusters 1 and 2 (Figures 4 and 5).Communities in cluster 4, representing severely polluted sites, had steeper RAD curves and a minimal number of species.All slope values ranged from ´0.12 (at YBV) to ´1.90 (at YPY).The absolute values of the slopes in clusters 1 and 2, representing the less polluted sites, fell within the range of 0.12-0.32,whereas the values in cluster 3 were relatively higher (within the range of 0.20-1.00).The slopes in cluster 4 were significantly higher than those of clusters 1 and 2, ranging from 0.64-1.90.

Indicator Species Analysis
Forty species were identified as indicators (p < 0.05) for four different clusters, according to the indicator species analysis (Figure 5).Among the indicator species, 22 were chosen in cluster 1 with an IndVal of 23.1-89.7,and 14 in cluster 2 with an IndVal of 34.5-71.5 (Figure 5).Meanwhile, three and one indicator species were identified in cluster 3 and cluster 4 (representing polluted streams), respectively.The Plecoptera species Nemoura KUa and Neoperla quadrata showed high indicator values for cluster 1. Ephemeroptera (Caenis nishinoae, Uracanthella rufa, and Ecdyonurus levis) and Tricoptera species (Hydroptila KUa, Cheumatopsyche KUa, and Hydropsyche kozhantschikovi) showed high IndVal (>50) in cluster 2. No aquatic insect was selected as indicator species in clusters 3 and 4. Tubificidae species were considered indicator species for cluster 4 that represented polluted sites (Figure 5).
According to WARMs, the indicator species in each cluster displayed different preferences toward land use types and three environmental variables (slope, conductivity, and substrate index) based on regression coefficients (Figure 5).Regression coefficients between the inferred and the observed proportions of land use types and environmental variables were relatively high, and were similar with and without bootstrapping tests (Table 3).This observation indicated that the optimal variables for each species were estimated reliably by the WARMs.The indicator species in cluster 1 showed a higher optimum preference in mountain forest areas with lower conductivity than in agricultural and urban areas with higher conductivity, ranging between 170-395 μs/cm.The indicator species, including Ecdyonurus kibunensis, Apatania KUa, and Tipula KUa, had high preferences in agricultural lowland areas with relatively high conductivity.Three species, Neoperla quadrata, Paraleptophlebia chocolata, and Baetis ursinus, had relatively high optimum preference in the urban areas of cluster 1.Although clusters 1 and 2 have similar RADs, their habitat preference scopes were substantially different (Figure 5).The proportion of forest area was generally lower in cluster 2 than in cluster 1, whereas that of agricultural area was higher in cluster 2. Species in cluster 1 preferred low conductivity and high slope, whereas the substrate preference was positive in clusters 1 and 2. Hydroptila KUa, Caenis nishinoae, Uracanthella rufa, and Eubrianax KUa showed high optimal preference values for the agriculture areas of cluster 2 (>40%).Hydropsyche kozhantschikovi and Hirudinidae species were well adapted to disturbed streams with high conductivity.In cluster 3, however, the preference was not much different from that of cluster 2, while the number of indicator species was substantially lower in cluster 3 than it was in cluster 2 (as stated above).In cluster 4, representing severe pollution, Oligochaeta species were exceptionally well adapted to urban areas with the highest conductivity and the lowest substrate size, showing a relatively high IndVal of 45.6 (Figure 5).

Indicator Species Analysis
Forty species were identified as indicators (p < 0.05) for four different clusters, according to the indicator species analysis (Figure 5).Among the indicator species, 22 were chosen in cluster 1 with an IndVal of 23.1-89.7,and 14 in cluster 2 with an IndVal of 34.5-71.5 (Figure 5).Meanwhile, three and one indicator species were identified in cluster 3 and cluster 4 (representing polluted streams), respectively.The Plecoptera species Nemoura KUa and Neoperla quadrata showed high indicator values for cluster 1. Ephemeroptera (Caenis nishinoae, Uracanthella rufa, and Ecdyonurus levis) and Tricoptera species (Hydroptila KUa, Cheumatopsyche KUa, and Hydropsyche kozhantschikovi) showed high IndVal (>50) in cluster 2. No aquatic insect was selected as indicator species in clusters 3 and 4. Tubificidae species were considered indicator species for cluster 4 that represented polluted sites (Figure 5).
According to WARMs, the indicator species in each cluster displayed different preferences toward land use types and three environmental variables (slope, conductivity, and substrate index) based on regression coefficients (Figure 5).Regression coefficients between the inferred and the observed proportions of land use types and environmental variables were relatively high, and were similar with and without bootstrapping tests (Table 3).This observation indicated that the optimal variables for each species were estimated reliably by the WARMs.The indicator species in cluster 1 showed a higher optimum preference in mountain forest areas with lower conductivity than in agricultural and urban areas with higher conductivity, ranging between 170-395 µs/cm.The indicator species, including Ecdyonurus kibunensis, Apatania KUa, and Tipula KUa, had high preferences in agricultural lowland areas with relatively high conductivity.Three species, Neoperla quadrata, Paraleptophlebia chocolata, and Baetis ursinus, had relatively high optimum preference in the urban areas of cluster 1.Although clusters 1 and 2 have similar RADs, their habitat preference scopes were substantially different (Figure 5).The proportion of forest area was generally lower in cluster 2 than in cluster 1, whereas that of agricultural area was higher in cluster 2. Species in cluster 1 preferred low conductivity and high slope, whereas the substrate preference was positive in clusters 1 and 2. Hydroptila KUa, Caenis nishinoae, Uracanthella rufa, and Eubrianax KUa showed high optimal preference values for the agriculture areas of cluster 2 (>40%).Hydropsyche kozhantschikovi and Hirudinidae species were well adapted to disturbed streams with high conductivity.In cluster 3, however, the preference was not much different from that of cluster 2, while the number of indicator species was substantially lower in cluster 3 than it was in cluster 2 (as stated above).In cluster 4, representing severe pollution, Oligochaeta species were exceptionally well adapted to urban areas with the highest conductivity and the lowest substrate size, showing a relatively high IndVal of 45.6 (Figure 5).

Community Patterns and Indicator Species
The applications of multivariate analysis, RADs, and indicator analysis effectively accounted for the environmental impact of human activity on benthic macroinvertebrates.Cluster analysis and NMS were useful to determine the changing patterns of benthic species composition influenced by physicochemical gradients and anthropogenic disturbance.Structural properties of communities described by RAD slopes, however, were feasible for assessing disturbance due to urbanization, and these results provided an additional dimension of biological responses in characterizing environmental impacts for monitoring ecosystems [21,26].Indicator species analysis using Indval and WARM determined the habitat preference and tolerance range across an environmental gradient [45].
Differentiation patterns on multivariate analyses can be further explained according to the properties of indicator species.Indicator species in each cluster presented higher values of occurrence rates and densities at sites in the same cluster than in other clusters.This indicated that indicator species with higher Indval play a key role in the community differentiation by multivariate analysis.Habitat preference range of indicator species based on WARM also provided additional information for patterning communities responding to different environmental impacts.For instance, the sampling sites belonging to clusters 1 and 2 were located on the left side of NMS axis 1, and had similar RAD patterns.Clusters 1 and 2, however, were clearly differentiated based upon habitat preferences.Species in cluster 2 had a relatively high preference for agricultural areas with low conductivity, whereas species in cluster 1 showed higher preferences for forest areas.It can be noted that community composition and structure are determined by indicator species analysis.
Indicator values provide an additional dimension for addressing ecological integrity in response to different environmental impacts.Community structure patterns based on RADs also show the relations between biological responses and environmental factors linked with NMS.Although numerous studies have reported species abundance distribution, few were concerned with integrating community structure, RADs, and indicative species.
Several species abundance distribution models have been proposed to identify patterns in community structures, including log-normal distribution, log series based on statistical models, and geometric series to show biological processes such as niche subdivision [27,49].Various parameters of these models (e.g., k values in the geometric series, a and γ values in the log-normal distribution) have been used as indicators for community structure [28,49,50].In this study, we extended the scope of RADs so that it would offer more information on community structure [21,27].The application of RAD slopes was supported by other studies that identified community responsiveness to environmental variability [26,41,42].This approach could also be used for the health assessment of various ecosystems.

Influence of Land Use Types
Patterns of land use in riparian areas can influence the suitability of habitats for stream communities, and change the distribution and abundance of species [51,52].In this study, most environmental variables, except depth, showed statistically significant correlations with forest and urban areas in the land use types (Table 1).Differences in land use types influenced biological properties, including community patterns (Figures 2 and 3 Table 2), indicator species (Figure 5) and environmental properties (Table 1).The proportions of urban and forest areas, which had high correlations with the first and second axes on NMS, and relatively strong correlations with RAD slopes, were considered as key factors in the maintenance of biodiversity and community stability overall [26].Similar observations were reported for numerous streams and rivers [53,54], confirming that the pattern of land use is a critical factor in the conservation of freshwater communities [11,55,56].
The land use type may affect biological trait patterns of benthic communities, such as functional feeding groups.Some shredders, such as Plecoptera and Trichoptera species in mountain streams, would be sensitive to organic pollution [57,58] and could serve as bioindicators in undisturbed forest areas, whereas filter feeders, like caddisflies and Tubificidae species (non-insect gatherers/collectors), may function as indicators for agricultural areas and urban streams contaminated with organic pollution, respectively [59].Further studies are warranted to learn more about the biological properties (e.g., traits) of benthic macroinvertebrate communities and their adaptation to land use types.

Geological Features and the Management of Urban Streams
Conserving biodiversity is possible if ecosystem management planning is conducted efficiently, even when the sample sites are located in the vicinity of human settlements.In this study, sampling sites such as YOU, YGS, YJJ, YSB, CCM, and CCD are located in the urban residential areas with apartment complexes and high population densities.These sites, however, displayed high biodiversity and biological indices, and were classified in cluster 1, which represents less polluted streams (Figure 2, Table 2).This phenomenon was the result of an active urban ecosystem management program initiated by the local government, such as the Stream Health Evaluation project, Dong-Cheon Stream Restoration project, or similar plans.These projects concentrated on the restoration and conservation of stream habitats resulting in increased habitat diversity, removing artificial structures, or reviving connectivity between streams and land [60].
However, some study sites were found to be polluted despite the fact that they were located in relatively less polluted areas.The study sites in clusters 2 and 3, which have relatively low biodiversity and biological indices, are in fact surrounded by mountains.This indicates that the areas were not properly ecologically managed.Specifically, these areas were disturbed by organic pollution due to excessive commercial development (DDK and DAU), sedimentation from cemetery construction (YPU), disturbance by local residents (JAU), and organic pollution due to agricultural land use (ZGU and TSD), etc.In the Busan region, human settlements occur mostly near the streams since the mountains are dispersed and there is little lowland area.Considering that there is a high risk of biological organisms being exposed to various anthropogenic pollution sources in aquatic ecosystems, these results demonstrate that urban ecosystem management plans should be carefully conducted to efficiently protect biodiversity and ecological integrity [56,61].

Conclusions
We investigated the changes on benthic macroinvertebrate communities in streams of urban areas.The changes of land use types in riparian zones due to the urbanization influenced on the species composition as well as community structure.RADs were effective for reflecting the overall degradation of community structure due to urbanization, and indicator species analysis provided additional information for identifying the gradient of habitat preference of species.These results suitably presented the ecological state of streams in urban areas in a comprehensive manner, including particular patterns such as the positive status of stream ecosystems in urban residential areas due to ecosystem management, or negative status in mountainous areas.These approaches would be suitable in assessing community structure properties and species composition, enabling a comprehensive evaluation of the ecological integrity of urban streams.

Figure 1 .
Figure 1.Location of sampling sites for each cluster with altitude (a); and land use type (b).

Figure 1 .
Figure 1.Location of sampling sites for each cluster with altitude (a); and land use type (b).

Water 2016, 8 , x 7 of 18 Figure 2 .
Figure2.Classification of sampling sites with densities of benthic macroinvertebrates through a hierarchical cluster analysis, using the Ward linkage method with Bray-Curtis distance.

Figure 2 .
Figure2.Classification of sampling sites with densities of benthic macroinvertebrates through a hierarchical cluster analysis, using the Ward linkage method with Bray-Curtis distance.

Figure 3 .
Figure 3. Nonmetric multidimensional scaling (NMS) ordination based on macroinvertebrate communities with fitted vectors of environmental variables (R 2 > 0.07).Different colors and symbols represent the four clusters.

Figure 3 .
Figure 3. Nonmetric multidimensional scaling (NMS) ordination based on macroinvertebrate communities with fitted vectors of environmental variables (R 2 > 0.07).Different colors and symbols represent the four clusters.

Figure 4 .
Figure 4. Rank abundance diagram (RAD) curves with mean and standard deviation (error bar) of each ranked species for three different clusters.

Figure 4 .Figure 5 .
Figure 4. Rank abundance diagram (RAD) curves with mean and standard deviation (error bar) of each ranked species for three different clusters.Water 2016, 8, x 11 of 18

Figure 5 .
Figure 5. Optimal and tolerance ranges of (a) land use types and (b) environmental variables for the selected indicator species.Dot indicates the optimal value, and the solid line indicates the tolerance value; * IndVal: indicator value, ** Occurring rate is frequency of occurrence (% of sites).

Table 2 .
Differences in parameters (mean ˘standard error) for different clusters defined by clustering analysis.
* p values indicate significant differences among clusters based on Kruskal-Wallis test.** The different alphabets indicate significant difference according to the Dunn's post hoc multiple comparison tests.

Table 3 .
Predictive power of weighted averaging regression models for land use types and environmental variables.
* R 2 : regression coefficient between the inferred and measured values.** RMSE: root mean squared error.

Table 3 .
Predictive power of weighted averaging regression models for land use types and environmental variables.
* R 2 : regression coefficient between the inferred and measured values.** RMSE: root mean squared error.