Spatial Distribution of Benthic Macroinvertebrate Assemblages in Relation to Environmental Variables in Korean Nationwide Streams

Conserving and enhancing freshwater biodiversity are global issues to ensure ecosystem integrity and sustainability. To meet this, it is critical to understand how the biological assemblages are determined by environmental gradients in different spatial scales. Nevertheless, information on their large-scale environmental relationships remains scarce in Korea. We aimed to understand nationwide spatial distribution patterns of benthic macroinvertebrates and important environmental factors affecting their distribution in 388 streams and rivers across Korea. A total of 340 taxa, belonging to 113 families in 23 orders of five phyla, were identified. Assemblage composition in most Korean streams included a few predominant colonizers and a majority of rare taxa. Cluster analysis based on benthic macroinvertebrates classified a total of 720 sampling sites into five clusters according to the pollution levels from fast-flowing less polluted streams with low electrical conductivity to moderately or severely polluted streams with high electrical conductivity and slow water velocity. Canonical correspondence analysis revealed that altitude, water velocity and streambed composition were the most important determinants, rather than watershed and water chemistry variables, for explaining the variation in macroinvertebrate assemblage patterns. The results provide basic information for establishing the conservation and restoration strategies of macroinvertebrate biodiversity against anthropogenic disturbances and developing more confident bio-assessment tools for diagnosing stream ecosystem integrity.


Introduction
Freshwater ecosystems occupying a very tiny fraction of the Earth's surface support remarkable biodiversity and abundance of benthic macroinvertebrates [1,2]. However, growing human pressures have substantially deteriorated the overall ecological integrity and induced a biodiversity crisis. This problem now increasingly becomes a global issue, and thereby calls attention to establish relevant conservation strategies and wise management practices to maintain sustainable freshwater environments.
Although Asian riverine ecosystems contain remarkable taxonomic diversity as well as high levels of endangerment and endemism, studies on biodiversity and ecology of these ecosystems have been poor [1,3]. Year-to-year variations in seasonality and its relationship to monsoonal climate, particularly in temperate Asian regions, can be major drivers that profoundly affect the hydrologic regime and geomorphology in stream environments, consequently determining the distribution and abundance of macroinvertebrates [4]. For example, a wet season low and a dry season high are expected for periodic seasonal patterns in abundance, depending on the frequency and intensity of summer monsoon rainfall. However, habitat destruction and water quality degradation in Asian river systems have become more dramatically epidemic than any other continents due to the fast growing human population and demands for economic development [5,6]. Because failure of conservation efforts has mainly resulted from scant evidence on river ecology, the priority over scientific and practical challenges to overcome this circumstance must be put on establishing species inventories and ecological information based on the accumulation of region-specific case studies, particularly at watershed or national scales [4,7].
Macroinvertebrate diversity and abundance are significant community attributes that are controlled by a variety of mechanisms at different spatial scales. A number of studies have documented how macroinvertebrate assemblages respond to environmental variables and which variables best explain their distribution and abundance. Some studies showed good relationships among macroinvertebrate assemblages, chemical variables [8], and the organic energy base [9], whereas habitat-related physical factors were widely demonstrated as primary contributors such as substrate composition [10], flow and current velocity [11], elevation and stream size [12] and temperature [13]. Vegetation, geology and human land use are also important for their spatial distribution [14,15]. While there remains a large number of documents focusing on such environmental relationships at finer spatial scales, studies on large-scale spatial patterns have been relatively small (e.g., at national scale [16,17] and at continental scale [18,19]).
Studies on the spatial distribution patterns of macroinvertebrate assemblages based on their environmental relationships are crucial. Firstly, those studies provide fundamental information for the conservation and restoration of biodiversity against anthropogenic disturbances. Secondly, knowledge of the species response to environmental gradients is important to separate the effects of pollution from the effects of natural variables on community structures [20]. Thirdly, it is possible to develop more confident bio-assessment tools for diagnosing stream ecosystem integrity because macroinvertebrate responds sensitively to environmental alterations that lead to changes in composition and community structures [21,22]. Thus, we performed a synoptic study of the large scale spatial distribution of benthic macroinvertebrates in relation to environmental variables. Specifically, the objectives of this study were: (i) to characterize the spatial distribution and assemblage structures of macroinvertebrates; (ii) to identify environmental distinction of Korean stream ecosystems based on their assemblages; (iii) to determine major environmental variables that affect their distribution; and (iv) to provide methodological considerations to improve biomonitoring programs, particularly focusing on sampling procedures. A relevant scientific basis would then be established with these data for developing sustainable management, conservation practices, and a reliable biomonitoring program. The results of this study could also contribute to a better understanding of the spatial distribution of freshwater macroinvertebrates in poorly explored Asian regions.

Study Area
Samples were collected at 720 sampling sites from 388 streams and rivers on a nationwide scale in South Korea. South Korea is located between 37˝00 1 N latitude and 127˝30 1 E longitude, encompassing the southern half of the Korean Peninsula with an entire area of approximately 100,033 km 2 . The annual precipitation is 1308 mm, but there is substantial variation among seasons [23]: severe flooding events during the summer monsoon period, and base flow or even drought conditions in the other seasons. The entire stream system throughout South Korea comprises five watersheds, including five large rivers, their tributaries, and other small independent streams ( Figure 1). Land use types in each watershed generally display well-managed forests upstream, and agricultural and urban development from middle to lower regions. Most Korean streams have suffered from a variety of human activities that alter the physicochemical stream environments, particularly caused by channel modification and eutrophication [24].
Among 720 sampling sites, the largest number of sites belonged to the Han River watershed (HRW) (n = 320), followed by the Nakdong River watershed (NRW) (n = 130), the Geum River watershed (GRW) (n = 130), the Youngsan River watershed (YRW) (n = 76), and the Seomjin River watershed (SRW) (n = 64). It was assumed that such a nationwide survey covered the majority of stream types in Korea to understand how lotic macroinvertebrates are spatially distributed in relation to environmental factors. Field sampling was conducted during spring (May 2009) under base flow conditions because the highest benthic macroinvertebrate diversity was expected at that time. More information of the five major river watersheds in Korea can be found in Hwang et al. [25]. Geographic locations of the study sites. Five macroinvertebrate-based site groups were classified as G1a (n = 126, black circles), G1b (n = 199, white circles), G2a (n = 106, white triangles), G2b (n = 118, gray triangles), and G2c (n = 171, black triangles) based on Sørenson distance measure cluster analysis. Dark lines indicate the five large rivers and the light lines display their tributaries and small independent streams.

Measurement of Environmental Variables
Environmental variables were measured both in the field and in the laboratory to define their effects on benthic macroinvertebrate assemblages. Three categories of environmental variables were considered such as watershed-related regional variables, physical in-stream properties, and water quality elements.
Altitude and land-use types were considered as regional variables for watershed characteristics. The altitude at each sampling site was obtained using digital elevation data (Openmate Inc., Seoul, Korea). The proportion of prevalent land use types was determined for each sampling site using a topographic map (1:50,000).
Physical in-stream properties were measured during field surveys: (i) stream width for the distance from bank to bank at a transect representative of the stream channel; (ii) wetted stream width using a range finder (model LRM-1500M, Newcon Optik Inc., Toronto, ON, Canada); (iii) water depth of the vertical distance from the water surface to stream bottom; (iv) current velocity at riffles or gliding runs using a current meter (3000-LX, Swoffer Instruments, Inc., Tukwila, WA, USA), or calculated by the Craig method [26,27]; (v) percent substrate composition visually estimated as fine (<2 mm) and coarse-sized particles (ě2 mm).
Water quality variables such as pH, water temperature, dissolved oxygen (DO), turbidity, and electrical conductivity (EC) were measured with a multi-probe portable meter (e.g., YSI 6920, YSI Inc., Yellow Springs, OH, USA or Horiba U-22XD, Kyoto, Japan) at the center of each sampling stretch. Two liters of water samples were collected in sterilized plastic bottles at each site, kept in a container with ice, and transported to the laboratory. Laboratory measurements were conducted to determine the biochemical oxygen demand (BOD), total nitrogen (TN), and total phosphorus (TP) following Standard methods [28].

Sampling of Benthic Macroinvertebrate
A Surber sampler (30 cmˆ30 cm, 1 mm mesh) was employed to collect benthic macroinvertebrates. Samplings were conducted at fast-flowing riffle or gliding run habitats (in the case when suitable riffles were not available) within 100 m. Quantitative samples were taken from three randomly selected riffles at each site, pooled together in a 500 mL plastic bottle with 80% ethanol, and labeled. The sampling device and procedures followed the guidelines of the National Aquatic Ecological Monitoring Program (NAEMP), Korea [27].
After field sampling, all organisms were hand-sorted from detritus and inorganic material, and stored in 70% ethanol. Subsampling was permitted only for dominant taxonomic groups (e.g., Oligochaeta, Ephemerellidae, Chironomidae, and Hydropsychidae) with large numbers of specimens available in each sample. Each macroinvertebrate specimen was identified to the lowest possible taxonomic level (usually genus or species). However, several taxa with limited systematic information (e.g., Coleoptera and Diptera) were identified only to the family level. All individuals were counted and converted to individuals/m 2 .

Data Analysis
Differences in environmental variables and macroinvertebrate assemblages were identified among the five major river watersheds by the Kruskal-Wallis test, followed by Dunn's nonparametric multiple comparison test if they are significantly different (p < 0.05). Spearman's rank correlation analysis was used to indicate the relationships between the most dominant taxa and environmental variables. A cluster analysis was conducted to classify the benthic macroinvertebrate communities using the flexible beta method (beta =´0.25) with the Sørenson distance measure [29]. Samples were classified into clusters based on similarities in the community composition. The multi-response permutation procedure (MRPP, [30]), which is a non-parametric method to test for differences in assemblage structure among a priori defined groups [31], was conducted to evaluate differences among clusters. Differences in environmental variables among clusters were evaluated using the Kruskal-Wallis test and Dunn's nonparametric multiple comparison test. The Kruskal-Wallis test, Dunn's test and Spearman correlation analysis were performed with STATISTICA software (StatSoft, Inc., version 7, Tulsa, OK, USA).
Indicator species analysis (IndVal, [32]) was performed to evaluate potential indicator species in each cluster defined in advance. The indicator value of a species is the product of its relative abundance and frequency (ˆ100), ranging from 0 (no indication) to 100 (perfect indication) [33]. A perfect indicator of a particular group should be faithful and exclusive to that group, never occurring in other groups [29]. A species in a cluster that had an indicator value greater than in any other cluster was defined as good indicator species for that cluster in this study. A Monte Carlo method was performed to test the significance of indicator values.
Canonical correspondence analysis (CCA) was used to relate macroinvertebrate assemblages to environmental variables and to identify which environmental variables could best differentiate among the clusters [34]. Monte-Carlo simulations were carried out to verify whether variables exerted a significant effect (p < 0.05) on macroinvertebrate distributions. For the cluster analysis and CCA, those taxa with less than 0.2% of total abundance were excluded to minimize the effects of rare taxa. Taxa abundance data were transformed to log (x + 1) in both analyses to down-weight the effects of dominant taxa. Square-root-transformation was used for the environmental parameters expressed in percentages (e.g., type of land use and substrate composition), and the other variables were transformed to log (x + 1) except pH. After these processes all environmental variables were rescaled in the range of 0 and 1 based on the minimum-maximum range normalization [25]. Spearman correlation coefficients between scores and environmental variables were calculated to assist interpretation of changes in community profile using STATISTICA software (StatSoft, Inc., version 7, Tulsa, OK, USA). Cluster analysis, IndVal, MRPP and CCA were conducted using PC-ORD software (ver. 4.25, MjM Software Design, Gleneden Beach, OR, USA) [34].

Environmental Characteristics
Despite a large variation, physicochemical factors were significantly different among watersheds (Table 1). On average, the altitude, percent forest, and water velocity of the HRW and SRW were the higher than the other watersheds. The SRW streambeds were very heterogeneous with the highest percentage of coarse-sized particles, whereas those of the YRW contained a large amount of fine sediment. High altitude, fast water velocity and substrate complexity in both the HRW and SRW indicated good in-stream habitat conditions and potential to support high biodiversity and abundance of macroinvertebrates (Table 1a). The other three watersheds (NRW, GRW and YRW) were characterized by low altitude and a high proportion of agricultural land use.
Water quality parameters varied more prominently than physical variables among river watersheds (Table 1b). The average concentrations of DO, BOD, TN, TP, and turbidity were generally lower in the SRW than those in the other watersheds. EC varied widely depending on the sampling location, with the highest values near estuaries. Nutrient-related factors were particularly high in the GRW, YRW, and HRW.

Benthic Macroinvertebrate Assemblages
A total of 340 taxa, belonging to 113 families in 23 orders of five phyla, were identified in the five major river watersheds during the survey. Most of these were aquatic insects (272 species) including 144 Ephemeroptera Plecoptera Trichoptera (EPT; 62 mayflies, 24 stoneflies and 58 caddisflies) taxa and 35 Dipteran species. Taxa richness from the 720 sampling sites ranged from 0 to 49 with a mean of 14.4 (˘9.1) with the highest diversity at HRW and the lowest in YRW (Table 2a). In contrast, homogeneous streambeds and nutrient enrichment attributed to relatively low Shannon diversity index and high dominance index in both the NRW and YRW. The proportion of EPT taxa richness ranged from 40.1% (GRW) to 54.8% (SRW). One-way ANOVA revealed that taxa abundance significantly differed among river watersheds at p < 0.01, whereas there were no great differences for the relative abundance of each taxonomic group among watersheds (Table 2b). Table 1. Descriptive statistics of (a) regional and physical instream variables and (b) chemical variables in the five major river watersheds: the Han River (HRW), Nakdong River (NRW), Geum River (GRW), Youngsan River (YRW) and Seomjin River Watershed (SRW). Top and bottom lines of each variable indicate the average with standard deviation (in parenthesis) and the range, respectively. Kruskal-Wallis test (K-W) was performed with the environmental variables to compare the differences of each variable among river watersheds. The same small letters indicate no significant difference based on Dunn's multiple comparison tests. Abbreviations: DO, dissolved oxygen; BOD, biochemical oxygen demand; EC, electrical conductivity; TN, total nitrogen; TP, total phosphorus.      Korean stream ecosystems were characterized by a few predominant colonizers and a majority of rare taxa in macroinvertebrate assemblage composition. The most abundant and widespread taxa were a worm (Limnodrilus gotoi Hatai) and midge larvae (Chironomidae spp.) with their relative abundance of approximately 50% of total density throughout the whole river watershed. Other dominant species were mayflies (i.e., Baetis fuscatus (Linnaeus), Epeorus pellucidus (Brodsky) and Uracanthella punctisetae (Matsumura)) and netspinning caddisflies (Cheumatopsyche brevilineata Iwata, Hydropsyche valvata Martynov and Hydropsyche kozhantschikovi Martynov) among the different stream and river systems, most of which were dominant in somewhat nutrient-rich habitats at middle or lower streams. However, over 50% (195 taxa) of the fauna was present with low occurring frequency (less than 2% of all sites) and 90% with low abundance (less than 0.2%).
The environmental relationships of the dominant taxa were stronger with the physical variables (i.e., altitude, water velocity and streambed conditions) than with the chemical variables, among which water velocity was the most significant parameter. Consequently, significant positive relationships existed in most dominant taxa with peak abundance at a moderate velocity of 50-100 cm/s, particularly for E. pellucidus (r = 0.410, p < 0.001) and C. brevilineata (r = 0.435, p < 0.001) (Figure 2). Baetiella tuberculata (Kazlauskas), U. punctisetae and H. valvata tended to occur in their highest densities at the fast velocity (120-140 cm/s). No clear tendency was observed for Chironomidae spp.

Macroinvertebrate-Based Site Classification
The cluster analysis, based on the similarity in the benthic macroinvertebrate composition, largely classified the 720 sampling sites into two clusters and subsequently sub-clustered them into five groups. As a result, 126 sampling sites were included in Group 1a; 199 sites in Groups 1b; 106 sites in Group 2a; 118 sites in Group 2b; and 171 sites in Group 2c (Figure 1). MRPP validated these five groups with significant differences (A = 0.484, p < 0.001). The differences in environmental variables among the clusters are shown in Figure 3. The most important indicator taxa for each cluster are shown with their indicator values in Table 3.  Table 3. Indicator values (%) for the most important species (p < 0.05) in each cluster group. Monte Carlo tests (999 permutations) were used to assess the significance of each species as an indicator for the respective group (G1a-G2c). In total 50 species whose contribution to total density was higher than 0.2% were analyzed. Less important species were not shown in this Each cluster was clearly differentiated according to the differences of instream physicochemical conditions and geographical location of sampling sites. The cluster analysis discriminated less polluted streams with low EC and fast flowing water (Group 1) from moderately or severely polluted steams with high EC and slow water velocity (Group 2). Group 1a (G1a), congregating in the HRW, consisted specifically of mountainous upper streams although a tenth of this group was scattered over the other watersheds except the YRW. This group consisted of oligotrophic streams with distinguishing features of the highest altitude, the lowest BOD and nutrient concentrations, and the fastest water velocity. Additionally, the catchment was predominantly comprised of forested area. The best indicator species for this group were characterized as high sensitivity against organic pollution (e.g., Rhyacophila nigrocephala Iwata, Epeorus nipponicus (Uéno), Drunella aculea (Allen) and Hydropsyche orientalis Martynov). Group 1b (G1b) contained the largest number of sampling sites and was widely distributed throughout all river watersheds but mostly encompassing the agricultural and forested catchment. G1a and G1b closely resembled chemical environments, but steam sites belonging to G1b displayed mesotrophic condition with slightly higher BOD and nutrients than those of G1a. There existed the largest number of indicator species in G1b, among which U. punctisetae showed the highest indicator value.
G2a mostly included sites located in middle reaches of large rivers and their tributaries, particularly in the NRW and SRW. Although Group 2a (G2a) was similar in overall environmental characteristics to G1b, G2a was characterized with slightly lower altitude, higher EC and a lower water velocity when compared with G1b. Sampling sites suffering from poor water quality with organic degradation and nutrient enrichment were confined to Group 2b with the highest BOD, TN, and TP concentrations. These sites were influenced by the highest degree of agriculture and/or urbanization. The significant indicator species of this group were Hirudo nipponia Whitman, Chironomini sp. and L. gotoi, indicating high trophy and saprobity [35]. Water quality conditions in Group 2c (G2c) were as poor as those of G2b with the highest EC. Sites in this group were characterized with the lowest water velocity and the largest proportion of fine particles because G2c was gathered by streams adjacent to estuaries, large rivers, and dammed streams. The only three weak indicator species appeared in G2c.

Environmental Variables Affecting Macroinvertebrate Distributions
A CCA was performed to understand how macroinvertebrates were distributed along environmental gradients (Figure 4). Total variability explained in the species data was 15.9% (Table 4). The eigenvalues of the first CCA axis (0.281) and the second CCA axis (0.101) were significant (p < 0.01; 99 Monte Carlo permutation test). All three CCA results were significant based on a Monte Carlo permutation test (p < 0.01).  Table 4 and the codes for each taxon are shown in Appendix Table A1. Table 4. Spearman's rank correlation coefficients and probabilities of environmental variables and the CCA axes (** p < 0.01, * p < 0.05) (a), and summary of CCA results (n = 720) (b). All axes were significant based on Monte Carlo permutation procedures. The five clusters were well separated in the CCA ordination plot (Figure 4a). Both the G1a and G1b groups with good physicochemical environments were positioned on the right side of the ordination plot, and G2b and G2c were on the opposite side. However, the G2a group straddled both sides. Eight environmental variables (i.e., BOD, TN, TP, altitude, water velocity, % forest, % fine particles, and % coarse particles) had significant correlations with the first axis, among which altitude was the most significant contributor (r = 0.793, p < 0.01), followed by % coarse particles (r = 0.747, p < 0.01) and % fine particles (r =´0.744, p < 0.01) ( Table 4). All physical factors except for % fine particles were in the opposite direction from the chemical variables (Figure 4c). This result suggests that the decline in altitude reflected deterioration in water quality, accompanied by increases in organic material and nutrients. The heterogeneity of the macroinvertebrate habitats also decreased with the longitudinal gradient toward downstream. The second axis was negatively related with stream order (r =´0.658, p < 0.01) and DO (r =´0.467, p < 0.01). High positive scores with the first axis were denoted for the rhithronic (i.e., pertaining to the headwaters) and intolerant species (e.g., Drunella aculea (Allen), E. nipponicus, Glossosoma sp., R. nigrocephala Iwata and H. orientalis), whereas negative scores were observed for the potamic (i.e., pertaining to rivers) and tolerant taxa (e.g., Micronecta sedula Horváth, Physa acuta Draparnaud, Labiobaetis atrebatinus (Eaton) and Asellus hilgendorfii Bovalius) (Figure 4b).

Discussion
The Asian monsoon region is a global biodiversity hotspot suffering from increasing anthropogenic disturbances, but aquatic biodiversity and ecosystem integrity remain poorly explored [3]. This is the same situation in Korea, but recent establishment of the National Aquatic Ecological Monitoring Program [36] opened a new era to assess ecosystem health and biodiversity in Korea. Our work presented in this paper takes advantage of the opportunity of such a nationwide scale of survey. Although the five river watersheds in Korea exhibited differences in environmental conditions and macroinvertebrate taxa abundance, overall taxonomic composition at each watershed displayed little difference (Table 2). Instead, a considerable spatial variation in benthic macroinvertebrate assemblages accounted for the combination of both ultimate (e.g., altitude and the degree of land use) and proximate factors (e.g., flow, stream bed substrate, BOD, EC, and nutrients). This finding was in line with previous studies of benthic diatoms in Korea [25] and macroinvertebrates in other Asian countries [12,37], suggesting the importance of various multi-scale factors in structuring macroinvertebrate assemblages [38]. Our results provide basal information for the sustainable management and conservation practices of stream ecosystems.

Macroinvertebrate Taxonomic Composition
The macroinvertebrate assemblages in Korean stream ecosystems generally bear resemblance to those in tropical and other temperate streams at higher taxonomic levels [35,39]. However, temperate Korean streams were relatively rich in Ephemeroptera, Diptera, and Trichoptera when compared with tropical streams for the higher biodiversity of Gastropoda, Decapoda and other insect orders such as Odonata and Hemiptera [39,40] (Table 2). The taxonomic composition also included a large number of rhithronic fauna that prefer stony substrates. For example, Baetidae, Heptageniidae, and Ephemerellidae mainly dominated the Ephemeroptera, and Rhyacophilidae and Hydropsychidae composed the Trichoptera in this study.
We confirmed a total of 340 macroinvertebrate taxa in this study. Of the macroinvertebrate taxa, chironomid midge larvae and a small minnow mayfly (B. fuscatus) were extensively encountered throughout Korean stream environments with the highest occurring frequency with 94% and 76% of total sampling sites, respectively. Temperate Korean streams were also characterized with a few predominant colonizers and a majority of rare taxa in macroinvertebrate taxonomic composition. Moreover, only six cosmopolitan taxa were noticeably comprised of >60% of all macroinvertebrate samples, which were L. gotoi, Chironomidae spp., U. punctisetae, C. brevilineata, H. valvata and H. kozhantschikovi. These dominant taxa were found to be significant indicators of mesotrophic or polytrophic streams in our study (Table 3), consistent with a result that most Korean streams and rivers were degraded in both chemical and biological status [36]. On the other hand, most Korean streams were occupied by a great number of rare macroinvertebrate species with a small distribution range and/or low abundance, as was also demonstrated in other studies [41,42].

Environmental Relationships with Macroinvertebrate Distribution
Benthic habitats are complex, and a variety of environmental variables acting at multiple spatial scales regulate the composition and distribution patterns of stream macroinvertebrate assemblages in an exclusive or synergistic fashion [38,43]. We revealed that the variables associated with altitude and in-stream habitats best accounted for the largest amount of variability in our macroinvertebrate data set, supporting the more important determinants of local environmental factors than broad or regional parameters [16,44]. The importance and role of local environmental variables have also been highlighted in other aquatic communities: aquatic macrophytes [45], freshwater phytoplankton [46], benthic diatoms [44], intertidal macroinvertebrates [47], and fish [48]. However, a comprehensive understanding of multispatial scales is needed because of significant correlations between macrohabitat and microhabitat characteristics, depending on the relative size of the area studied [3].
In our study, the variability among the macroinvertebrate-based stream groups was more prominently explained by the altitudinal gradients together with streambed composition and water velocity than chemical variables (Table 3). First, the most widely accepted theory related to altitudinal changes is the river continuum concept (RCC), which displays structural and functional responses to the longitudinal gradient [49]. Altitude has also been well documented in other studies as a main descriptor to determine macroinvertebrate richness as well as other environmental variables such as temperature, hydrology, food availability, streambed condition, and water chemistry [3,50]. However, we revealed poor relationships of altitude with water chemistry in contrast with significant associations with physical variables and biological attributes such as taxa richness, EPT richness, EPT abundance, and the Shannon diversity index. This may be due to the fact that over half of the studied streams corresponded to lowland streams below 100 m a.s.l., which were characterized by moderate or slightly poor water quality with severe variations in BOD, TN, TP, and Shannon diversity-based saprobity [36,51] (Table 1). Harding et al. [52] demonstrated that increasing human land use intensity along a river continuum caused water quality degradation, consequently leading to changes in taxonomic composition from intolerant EPT dominated to tolerant taxa-dominated assemblages despite a gradual increase in taxa abundance by a few dominant species.
Riffle habitats are commonly characterized by shallow water depth, oxygenating fast-flowing water, and stony beds, and exhibit higher taxa richness and abundance than that at pools or habitats with fine sediment, as in our case [3,42,53]. Such hydraulic conditions are critical determinants for the distribution and species composition of benthic organisms [16,54] and are the driving forces for evolution of their morphologies and life history [11]. Rheophilous (i.e., having an affinity for running waters) and limnophilous taxa were clearly discriminated based on the correlation analysis ( Figure 2) and CCA (Figure 4) results. For example, well-known rhithronic species in two Ephemeropteran (Ephemerellidae and Heptageniidae) and one Trichopteran (Hydropsychidae) families mostly displayed high preference to the mid-or fast current conditions, which corresponds to the intolerant scraper or filtering collector groups abundant in upper and middle stream reaches [49]. Thus, biological monitoring and assessment programs using rheophilic macroinvertebrates would benefit from their ecological characteristics and their indication of good environmental quality (e.g., [42]) although seasonal fluctuations in hydraulic parameters by the Asian monsoon more critically cause a catastrophic drift and washout of benthic organisms.
Streambed composition is one of the most important factors to directly influence richness and abundance of macroinvertebrates on local scales [12,16]; thus, close correlations would be expected with variables related to longitudinal changes in stream ecosystems. Boulders and cobbles are typically the major structural elements in steep gradient upper streams, whereas sand and smaller sediments predominate in the lower reaches [55]. There have been a large number of studies on the macroinvertebrate-substrate relationship, most of which have revealed that macroinvertebrate diversity and density increase with higher heterogeneity due to the available stable and diverse microhabitats (e.g., [10,43]). Therefore, factors determining macroinvertebrate communities at local scales obviously put a priority on the streambed conditions rather than water quality or other physical variables. We found that the mainstreams of the Nakdong River, which contained about 60% fine particles, retained the lowest taxa richness (8.2 on average, n = 18) and abundance (747 individuals/m 2 ) among five major rivers despite its good water quality condition (mean BOD, 1.9 and TN, 1.8 mg/L) with the exception of the severely polluted Youngsan River (8.9 taxa richness and 7.3 BOD, n = 13). These results indicate that homogeneous streambeds with greater fine particles support lower diversity and abundance even when streams maintain good water quality [56]. Additionally, organic pollution could transform the coarse substrate into an organic rich soft bottom, altering the community structure from dominated by diverse and intolerant species to communities predominated by a few tolerant species [3,20], as in the case of the Youngsan River.

Considerations to Improve Macroinvertebrate Biomonitoring Programs
Not until 2006 did Korea adopt biological water quality criteria and the concept of ecological integrity in the water quality program [25,36]. Since then, Korean government has led a biological survey of streams and rivers (i.e., NAEMP) every year to assess the current biological status of stream and river ecosystems, and to develop a strategy for the restoration and management of disturbed systems [25]. As a part of a nationwide survey benthic macroinvertebrates are also monitored based on the guidelines and assessment tools [27]. Notwithstanding their suitability for convenient and rapid bio-assessment, there remain debatable issues as to whether the methods for sampling and treating macroinvertebrates are effective to provide a reliable and accurate indication of the macroinvertebrate fauna throughout the country.
The Korean national biomonitoring program presents field sampling in a cost-effective and time-saving manner for rapid bio-assessment. To satisfy this purpose benthic macroinvertebrates are optimized to be collected at riffle and/or gliding run habitats using a Surber sampler with three replicates [27]. However, such sampling methods probably underestimate overall biodiversity in stream ecosystems, considering the whole stream environment. First, single-habitat sampling possibly produces incomplete taxa lists and includes no target habitats at certain sites despite its advantage that the influences of both water and habitat quality on macroinvertebrates are not confounded by instream habitat variation [57]. In this regard, a multihabitat approach would be more profitable for the estimation of taxonomic diversity due to its consistent application across stream types especially at large-scaled survey, comprehensive taxa lists, and effective assessment of ecological conditions [58]. Second, the Surber sampler is one of the most commonly used quantitative tools in lotic systems and provides high-precision information on the abundance and composition of macroinvertebrate assemblages [21,59]. This method is, on the other hand, usually more appropriate for riffle habitats of shallow streams, presumably underestimating overall biodiversity in a region. Recent studies on the comparison between sampling devices suggested that artificial substrates (e.g., leaf-bags) would be used as complementary tools due to their discriminative taxonomic composition [60][61][62]. Third, a majority of biomonitoring programs widely adopt three to five replicates for lotic ecosystems because one way to reduce monitoring costs is to decrease the sample size [21,63,64]. The previous studies also showed that such a small number of replicates rarely influenced ecological health assessment, particularly for the indices applying sensitivity/tolerance taxa, as the same for NAEMP [65][66][67]. Nevertheless, small sample sizes may influence the values of biological measures and the representativeness for a real benthic community based on the asymptotic relationship of the number of taxa with both the sampling area and the collected individuals [68]. Finally, sieve mesh size could also affect the accurate estimates of taxonomic diversity and environmental quality assessment [67]. Finer mesh sizes more accurately represent macroinvertebrate assemblages than coarser mesh sizes, whereas they need more efforts to handle specimens. On the contrary, coarser mesh sizes undervalue the original density and taxa richness of macroinvertebrate assemblages by passing small individuals through the sieves and result in higher diversity [69]. Considering the tradeoff between cost and effectiveness, many rapid bio-assessments determine 0.5 mm mesh size as a reasonable choice [21,64,70]. Therefore, future researches are required for identifying optimal sampling effort comprehensively considering cost-effectiveness, easy applicability, and well representative of resident macroinvertebrate assemblages.

Conclusions
Large-scale knowledge on environmental relationships of aquatic communities is crucial to conserve freshwater biodiversity and sustain ecological integrity. Korean stream macroinvertebrate assemblages were determined not only by regional and physical instream variables but also by pollution-related parameters at nation-wide scale. Macroinvertebrate-based site classification evidently provided an environmental characterization for different types of Korean streams, indicating that benthic macroinvertebrates are a valuable biomonitoring material. The results of this study provide important information and a bridge for further work such as the assemblage-specific responses to environmental disturbance. Our results also contribute to establishing effective management practices, implementing conservation measures, and developing more reliable biological monitoring tools for sustainable freshwater ecosystems.
Author Contributions: Yung-Chul Jun developed the concept of the study, performed data analysis, and wrote the manuscript. Nan-Young Kim assisted data processing and figure construction. Sang-Hun Kim conducted quality assurance of the physico-chemical and biological data. Young-Seuk Park and Dong-Soo Kong assisted study design and statistical analysis. Soon-Jin Hwang supervised the research, and assisted data interpretation and manuscript preparation. All the authors contributed to the review of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Codes for macroinvertebrates contributing to more than 0.2% of total abundance for canonical correspondence analysis in Figure 4.