Optimization of the Korean river macroinvertebrate prediction and assessment System (KRIMPAS) for the biological assessment of rivers

Abstract The biotic community predictive model proposed in the River Invertebrate Prediction And Classification System which compares observed values at test sites with expected values in reference sites is used to perform biological assessments of rivers. In 2021, a Korean River Macroinvertebrate Prediction and Assessment System (KRIMPAS) for river integrity assessment was developed. However, since it was the early stage of development, improvement of the accuracy and precision of assessment was required. Therefore, benthic macroinvertebrates and six environmental variables collected from 357 reference sites between 2010 to 2021 were used to compare the prediction mothods and attempt to optimize the KRIMPAS. Based on the number of taxa (species and family levels) and six biological indices, the existing KRIMPAS method (Original method) was compared with two weighted average methods (W-Avg 1 and 2) according to the probability that a test site was assigned. In terms of statistical significance, observed to expected (O/E) values of the number of taxa, but, there were difference for most of biological indices. When the probability density function based on the log-normal distribution for each derivation method was applied, no difference was observed in the number of taxa; however, there were differences in the mean, median, mode, and variance for each biological index. A fit verification based on the normalized root-mean-square error (NRMSE) for each derivation method was conducted, and the overall fit tended to be higher with the W-Avg 1 and 2 when compared with the Original method. Thus, to correct the current biological assessments of Korean rivers we propose the application of the W-Avg method as this would improve the derivation of the predictive model.


Introduction
The methods used to assess the integrity of river ecosystems worldwide were recently divided into two main categories (Herman and Nejadhashemi 2015).The first category involves multimetric index (Herlihy et al. 2005;Meier 2006;Touron-Poncet et al. 2014; USEPA 2020) while the second includes O/E (expected to observed) systems (Davies, 2000;Yuan et al. 2008;USEPA 2013;Clarke and Davy-Bowker 2018).Both methods have commonalities, as they use the environmental variables and biotic community data of reference streams (or sites) with minimal anthropogenic stress (Poquet et al. 2009;USEPA 2013;Nichols et al. 2014;Silva et al. 2017).Multimetric indexes, however, use multiple biological indices representing various environmental variables to assess the integrity of streams, whereas O/E systems use an empirical model to assess sampling sites (or test sites) and compare the predicted and observed biotic communities.O/E systems have been applied for the biological assessment and monitoring of river ecosystems in the United Kingdom (UK), Australia, and other European countries (Clarke et al. 2003;Parsons et al. 2004;Poquet et al. 2009;Pardo et al. 2014;Consoli et al. 2021).While the United States of America uses a multimetric index as the main assessment method nationally, they have also developed and introduced an O/E system for the assessment of specific states or rivers (Hawkins et al. 2000;Carlisle and Meador 2007;Hargett et al. 2007;Yuan et al. 2008).
The UK's River Invertebrate Prediction and Classification System (RIVPACS) was the first O/E system developed based on the benthic macroinvertebrate community of rivers in each country or region (Moss et al. 1987).The assessments compare the expected biotic communities with the observed biological communities in the absence of artificial stress for aquatic environment, and the prediction results differ depending on the observed biotic community and environmental conditions measured at the reference sites (Clarke et al. 2011;Clarke and Davy-Bowker 2018).The O/E system is calculated as the ratio of the actual value to the predicted value (or expected value) at the sampling sites in the reference streams; therefore, the lower the O/E value, the more stressed the river (Aroviita;2009;Hawkins 2009;Yuan 2006).As the O/E ratio predicts the composition of the taxa that can inhabit a sampling site, it can also be used to make predictions regarding taxa known to be sensitive to stress, which can lead to an inferred diagnosis of the stressors most likely to affect the region (Chessman 2021;Consoli et al. 2021).In addition, O/E systems can also use separate biological indices to measure the stress of a specific environmental variable, regardless of the number of taxa, as well as a predictive model (Davy-Bowker et al. 2010;Clarke et al. 2011).This makes it possible to comprehensively assess the state of a river using indicators that represent the complexity of various environmental variables.
In the Republic of Korea, several biological indices based on benthic macroinvertebrates have been developed for the biological assessment and monitoring of rivers (Kong et al. 2013;Kong and Kim 2016;Kong, Park et al. 2018;Kong, Son et al. 2018;2019), including the Benthic Macroinvertebrate Index (BMI; Kong, Son, et al. 2018), which is currently used in the Nationwide Aquatic Ecological Monitoring Program (NAEMP) by Ministry of Environment (ME) and National Institute of Environmental Research (NIER; NIER 2019).However, as this index is based on 5-day biochemical oxygen demand (BOD 5 ), it can only assess the impact on a specific environmental variable, such as organic pollution, and the conditions of all rivers are assessed on the same basis without considering stream type (Kong and Kim 2016;NIER 2019).Min and Kong (2021) have applied the RIVPACS prediction methods and attempted to develop an Korean River Macroinvertebrate Prediction and Assessment System (KRIMPAS) which is suitable O/E system for river ecosystems in the Republic of Korea.Their study first selected reference sites by considering stream types and proposed a prediction methods for each species and family level that was appropriate for the situation at the time of the early model development.
The biotic communities predicted by RIVPACS are determined probabilistically by suggesting taxa that can inhabit an area when there is no or minimal stress, based on the measured environmental variables at the sampling site, and the index is based on the predicted taxa (Davy-Bowker et al. 2010;Clarke et al. 2011;Clarke and Davy-Bowker 2018).The prediction of the probability of occurrence or individual abundances for each taxon is reflected as a weighted average value using the probability of belonging to each reference group based on the frequency or mean individual abundance of each taxon within the group, and the measurement values of the environmental variables of the sampling site (Moss et al. 1987;Yuan 2006;Clarke et al. 2011).If it is assumed that there are innumerable sampling units available, the distribution of the probability of occurrence and individual abundance for each taxon predicted based on environmental variables assumes asymptotic normality according to the central limit theorem (Goldstein 2009;Bell 2015).This eventually shows the distribution of the probability of occurrence or individual abundance with high frequency based on the mean value of each taxon.Therefore, the differences in the predicted index values for each sampling site may not be large, and the consequences of incorrect predictions can lead to errors in assessment (Clarke and Davy-Bowker 2018;Chessman 2021).For example, if the predicted value of the index is low, it can lead to an overestimation for the upper stream area with less stress; conversely, if the predicted value is high, the conditions downstream can be underestimated, even if it is a reference river.The purpose of this study was thus to predict the number of taxa and values of each index using the mean values for each group at different reference sites and the probability that sampling sites are assigned to each reference group, as well as to compare them with existing methods (deriving values of each index based on predicted taxa).The resulting data should also help facilitate the proposal of a more appropriate prediction methods for current use.

Study area
This study has involved the five major rivers distributed in the Republic of Korea, namely the Han, Nakdong, Geum, Yeongsan, and Seomjin watersheds.Data on benthic macroinvertebrates and environmental variables was thus collected from 357 reference sites (2,672 sampling units) located in the main streams, tributaries, and other independent streams of these rivers (Figure 1).Approximately 60% of the Republic of Korea consists of mountainous terrain (Kim et al. 2017).In the eastern region there is a concentration of high-altitude mountains which form a steep slope.The altitude decreases in the western region which has a high proportion of low-lying plains (Lee et al. 2011;Li et al. 2012).Based on these topographical characteristics, mountain streams at high and low altitudes were found to be more readily formed than plain rivers (Bae et al. 2003).Additional information on the rivers in the Republic of Korea was provided by Hwang et al. (2011), Jun et al. (2012), and Min and Kong (2021).

Data collection
This study used the results of the NAEMP conducted by the ME and NIER from 2010 to 2021.Environmental variables and benthic macroinvertebrates were sampled according to NAEMP guidelines (NIER 2019).
Quantitative sampling of benthic macroinvertebrates was conducted using three replicates per site (0.27 m 2 for each sampling site) using a Surber sampler (30 cm × 30 cm, 1 mm mesh; NIER 2019).Sampling was conducted annually in spring and autumn (some sites were surveyed once a year owing to inaccessibility or dry streams), focusing on the riffles at each site, and as stream size increased and the riffles became absent or difficult to access, the runs and pools were targeted instead.The benthic macroinvertebrates collected at each site were placed in a collection bottle, fixed in 95% ethanol, transported to the laboratory, and stored in 80% ethanol until further sorting and identification.All specimens were identified at the species level where possible; however, taxa for which taxonomic information is limited (e.g.Oligochaeta, Hemiptera, Coleoptera, and Diptera) were classified at the genus or family level.The individual abundance of each taxon for which the sample analysis was completed was converted to individual abundance/m 2 .
A total of nine environmental variables (catchment area, altitude, latitude, longitude, water depth, velocity, composition of the streambed, BOD 5 , and TP) were used in this study, and they were measured prior to the benthic macroinvertebrate sampling at each site.Environmental variables were used to confirm the environmental characteristics of each reference group proposed by Min and Kong (2021), and to derive the predicted values for the number of taxa and biological indices using a predictive model (KRIMPAS).The catchment area representing the stream size was derived from data from the Korean Reach File (KRF) provided by NIER (2015).Altitude, latitude, and longitude data representing the geographical characteristics were measured using a GPS device (Triton 500, Magellen Inc., San Dimas, CA, USA).For the water depth, velocity, and composition of the streambed representing the microhabitat, the values at the sampling points where the benthic macroinvertebrates were collected were measured and then averaged.Water depth was measured as the vertical distance from the water surface to the streambed.Velocity was measured using a flowmeter (Flow-Mate Model 2000 or 3000-LX, Swoffer Instruments, Ins., Tukwila, WA, USA) or based on the method proposed by Craig (1987) (1962), and the area occupied by the ratio of each substrate in the points was visually measured (NIER 2019).Subsequently, the mean diameter of the substrate of each streambed was derived from the weighted average of the median particle size of each substrate (boulder, −9 ɸ; cobble, −7 ɸ; pebble, −5 ɸ; gravel, −2.5 ɸ; and sand, 0 ɸ) and the area ratio occupied by the sites.The 5-day biochemical oxygen demand (BOD 5 ) and total phosphate (TP) corresponding to the water quality were derived from the water sampled at the sites according to the water pollution process test standards (NIER, 2021).

Reference sites, river groups, and predictive models
This study aims to compare and review the prediction methods of previously developed predictive models.Therefore, the 357 reference sites in the Republic of Korea (Figure 2 and 3) proposed by Min and Kong (2021) and the three predictive equations used to derive the predictive model were used uniformly.In Min and Kong (2021), reference sites were selected using six stream types and five environmental variables for 2,700 sites located in the Republic of Korea (Table 1).The selected 357 reference sites were classified into six reference groups by Two-Way Indicator Species Analysis (TWINSPAN) based on benthic macroinvertebrate communities.Stepwise discriminant analysis (MDA) was performed to extract environmental variables that differentiate each reference group.As a result, a total of five discriminant functions and six environmental variables (catchment area, latitude, velocity, water depth, altitude, longitude) and a predictive model was derived using discriminant functions 1 to 3 (predictive equations) among them.At this time, as available data were added and compared with those of Min and Kong (2021), environmental variables and benthic macroinvertebrates by site were newly obtained and applied as mean values for the data used in this study.
Groups 5 and 1 and the highest and lowest number of sites in their reference groups, respectively (Figure 2).From the perspective of each watershed, the sites in Groups 1 and 2 were mainly distributed in the Han and Nakdong Rivers, while most of the sites in Group 3 were in the Han River.The sites for Group 4 were evenly distributed throughout the four watersheds, except for the Yeongsan River, whereas the sites of Group 5 were mainly located in the Nakdong and Seomjin Rivers.Most sites in Group 6 were located on the Nakdong River.
Based on the mean value of each environmental variable, the stream size was found to be small in Groups 1, 2, and 5, and largest in Group 6 (Figure 3).Groups 1 and 6 had the highest and lowest altitudes, respectively, and this variable tended to decrease overall in the order of Groups 1-6.The latitude and longitude changed depending on the position of the sites, and were highest in Group 3, where most of the sites were in the Han River, and lowest in Groups 4 and 5, where sites were evenly distributed over the entire watershed or concentrated in the Nakdong and Seomjin Rivers.There was little difference in water depth between Groups 1 and 5, but it was relatively deep in Group 6 when compared to the other groups.The velocity showed a decreasing trend in the order of Groups 1 to 6, like the altitude, and the decrease was most rapid in Group 1 and slowest in Group 6.The mean diameter of the streambed also showed similar changes with altitude and velocity, with a high proportion of coarse particles in Groups 1-6.Similar to the catchment area, TP tended to be low in Groups 1, 2, and 5, and high in Groups 3, 4 and 6.BOD 5 was low in Groups 1, 3, and 4, and high in Groups 2, 5, and 6, so did not show similar trends to other environmental variables.

Data analysis
Three methods were reviewed in this study.The first was the original KRIMPAS-based method (Ori) proposed by Min and Kong (2021) which only used the BMI indicator taxa at the species level and the Benthic Macroinvertebrates family-level biotic Index (BMFI; Kong et al. 2019) indicator taxa at the family level to predict, observe, and derive O/E values.The second was the mean values (number of taxa and indices) for each reference group which were used to derive the predicted values using the weighted average method (W-Avg 1) based on the probability that a site was assigned to each reference group, in which the entire taxon list (all taxa) was used to derive predicted and observed values.The third method was the same as the second method, except it only used the indicator taxa of BMI or BMFI (W-Avg 2), as in the first method (Ori).Therefore, the first and third methods calculated O/E based on predicted and observed values using only indicator taxa, whereas the second method derived predicted and observed O/E values using all taxa.
When a total of eight metrics were used, three predicted values (Ori, W-Avg 1 and W-Avg 2) and two observed values (all taxa and indicator taxa) were derived for the number of taxa (species and family levels, Table 2).For biological indices, two predicted values (Ori and W-Avg 1) and one observed value (all taxa) were derived.The number of taxa (species and family level)  and the Total Ecological Score of Benthic macroinvertebrate communities (TESB; Kong, Park, et al. 2018) confirmed the community diversity within the sites.BMI, Average Ecological Score of Benthic macroinvertebrate community (AESB; Kong, Park, et al. 2018), and BMFI are the statuses for biological water quality based on BOD 5 .The Benthic Macroinvertebrates Streambed Index (BMSI; Kim and Kong 2019) is the status of the physical environment based on the composition of the streambed, and the Korean Thermality Index (KTI; Kong et al. 2013) is an index developed to assess the distribution of cold-water organisms in response to climate change.At this time, all three methods of prediction and two methods of observation can be applied to the number of taxa; however, when only the indicator taxa of BMI and BMFI are used in BMSI, KTI, TESB, and AESB, incorrect index values may be derived because of the omission of the indicator taxa required for each index (Appendix 1-5).Therefore, the predictions for the six indices used all taxa, but the predicted values were derived by applying Ori and W-Avg 1, and only the results derived when using all taxa were used for the observed values for comparison.
To compare the distributions of the O/E derived from the predicted and observed values according to each method, the probability density function and cumulative distribution functions for the lognormal distribution model were used based on the data from the reference sites.To derive additional statistical differences, one-way analysis of variance (ANOVA) was performed for the number of taxa (species and family levels) and independent-samples t-tests was conducted for the biological indices.The one-way ANOVA was performed for the number of taxa, and the t-tests were performed for BMI, BMSI, KTI, TESB, AESB, and BMFI.
The Pearson correlation analysis between the Ori, W-Avg 1, and W-Avg 2, and the degree of consistency between the newly proposed method (W-Avg 1 and 2) and the existing method (Ori) were confirmed.Finally, the fit of the predictive model for each method was reviewed based on the normalized root mean square error (NRMSE) of the predicted and observed values.
One-way ANOVA, t-tests, and Pearson's correlation analyses were performed using PASW (for Windows 8, version 6.3; SPSS Inc., Chicago, IL, USA).Variables for which the raw data did not show a form close to normality were converted into ln( ) x +1 .

Predicted, observed, and O/E values for each reference group
By applying Ori, W-Avg 1, and W-Avg 2 to each reference group, the predicted and observed values for each metric were obtained, and the O/E derived (Figures 4 and 5).The predicted and observed values for the number of species or families corresponding to the number of taxa tended to be higher with the W-Avg than the Ori for all reference groups, but there were no differences in the methods used for O/E values.
In the process of deriving the six metrics corresponding to the biological indices, the indicator taxa used for each index were slightly different, so it was impossible to confirm the differences in the observed values using the taxa methods (using indicator taxa or all taxa).The BMI, BMSI, AESB, and BMFI tended to have higher predicted values for W-Avg 1 than for Ori in the overall reference group.While the difference in the results were not notably large for BMI and AESB, there were cases in which the assessment grades were different.For BMSI and BMFI, the predicted values for the W-Avg 1 were higher than those for the Ori for all reference groups, and the differences in values were larger than those for the BMI and AESB, and consequently, there was a difference in the assessment grade.In the O/E for the BMI, BMSI, AESB, and BMFI, the Ori's results tended to be relatively high in contrast to the predicted values, but the W-Avg 1 values were generally close to 1.0 in most reference groups.Kong, son, et al. 2018) Benthic Macroinvertebrates streambed Index (BMsI; Kong and Kim 2016;Kim and Kong 2019) Korean thermality Index (KtI; Kong et al. 2013) total ecological score of Benthic macroinvertebrate community (tesB; Kong, Park, et al. 2018) average ecological score of Benthic macroinvertebrate community (aesB; Kong, Park, et al. 2018) Benthic Macroinvertebrates family-level Biotic Index (BMfI; Kong et al. 2019) For the KTI and TESB, the results from the Ori were higher when compared with the W-Avg 1 for all reference groups.The KTI did not show a large numerical difference between the Ori and W-Avg 1, but there were cases in which there were differences in grades according to the assessment results.For the TESB, the numerical differences between the derivation methods was very large; however, due to the nature of being affected by the number of taxa that appeared, in terms of assessment grade, there were no differences with all A's.In the O/E, the W-Avg 1 showed high values than the Ori for all reference groups, and when the O/E for each group was close to 1.0 so too was the W-Avg 1.The O/E for the Ori of the TESB was < 0.2, and thus very different from 1.0.

O/E Differences between the original and new methods for the reference sites
The distributions of O/E for Ori, W-Avg 1, and W-Avg 2 for the 357 reference sites were confirmed using a lognormal distribution (Table 3 and Figure 6).There was little difference in the mean, median, mode, and variance between the three methods for the number of species and families.Most of the sites ranged from 0.8-1.2 for their O/E values and showed a form of positive skewing.Among the metircs applicable to biological assessments, the BMI and AESB showed little difference in the mean, median, and mode between the Ori and W-Avg 1, similar to the number of taxa, but there was a slight difference in variance due to the difference in the number of sites corresponding to 0.9-1.2O/E values.
For the BMSI and BMFI, the mean, median, mode, and variance tended to be higher in the Ori than in W-Avg 1.In addition, in the BMSI, most of the sites were distributed in the O/E range of 1.1-1.5 for the Ori, whereas for the W-Avg 1, sites were distributed intensively in the 0.9-1.2range, so the mode W-Avg 1 mode was relatively close to 1.0.BMFI also showed a similar distribution to BMSI, and most of the sites were concentrated in the O/E range of 1.1-1.4 for the Ori and 0.9-1.2 for the W-Avg 1.The two indices had positively skewed distributions close to the normal distribution for all methods.
In the KTI and TESB, the mean, median, and mode tended to be higher in the W-Avg 1 than in the Ori, indicating opposite tendencies, and in terms of variance, there was little difference in KTI, but a large numerical difference was observed in TESB.In the KTI, sites were mostly distributed in O/E range of 0.7-1.2 for the Ori, whereas for the W-Avg 1, they were concentrated in the O/E range of 0.8-1.3.Among them, as most sites ranged from 0.9-1.0, the mode for the W-Avg 1 was relatively close to 1.0.For the TESB, sites were mostly distributed in the O/E range of 0.1-0.2 for the Ori, whereas for the W-Avg 1 the O/E ranged from 0.6-1.4.Therefore, the difference in the mean, median, mode, and variance between the Ori and W-Avg 1 were very large when compared to the other indices, and the mode for the W-Avg 1 was relatively close to 1.0.All methods showed a positively skewed form.
Based on the data from the reference sites, one-way ANOVA was applied to the number of taxa from which the the O/E values for the three methods (Ori, W-Avg 1, and 2) could be derived, there was no significant difference at either the species or family levels (Table 4).For BMI, BMSI, KTI, TESB, AESB, and BMFI, from which O/E for the two methods (Ori, W-Avg 1) could be derived, the O/E difference between each method was confirmed using a t-test.(Table 5).There were no statistically significant differences in BMI and AESB; however, significant differences were observed at a levels of < 0.001 for the BMSI, KTI, TESB, and BMFI.

O/E Correlation and fit between original and new methods
Pearson's correlation analysis was performed to confirm the correlation between the O/E for the Ori and W-Avg (1 and 2) for the 357 reference sites (Table 6).The correlation coefficients between the same metrics were all > 0.8, and the p-value were all < 0.01, indicating significant results.There was no significant difference in the number of taxa for the O/E between the three methods, but the Ori had a relatively high correlation with W-Avg 2, which used only the indicator taxa for BMI and BMFI, when compared to W-Avg 1, which used all taxa and showed a correlation coefficient of 1.0.In addition to the same metric, a strong correlation was observed between metrics with similar characteristics.The number of species and families and TESB, which were most affected by the presence or absence of taxa, showed higher correlation coefficients than the other metrics.BMI, AESB, and BMFI, which represent the degree of organic pollution, also showed relatively higher correlation coefficients than the other indices.In when the p-value of leven's equality of variance test was higher than 0.05 and the tukey hsD was lower than 0.05, the games-howell post hoc test was applied for the one-way anoVa p-value addition, BMSI, which represents the composition of the streambed, and KTI, which represents the distribution of the cold-water taxa based on altitude, showed relatively high correlation coefficients with BMI, AESB, and BMFI, even though they were developed based on environmental variables unrelated to organic pollution.Assessments of the correlation between the number of taxa and TESB, showed that the correlation between Ori and W-Avg 1 and 2 tended to be relatively higher than that between W-Avg 1 and 2 and the Ori methods.For example, the number of Ori species showed a correlation coefficient of 0.983 with W-Avg 1 TESB, whereas the number of species for W-Avg 1 and 2 showed correlation coefficients of 0.819 and 0.818 with the Ori's TESB, respectively, and the former was relatively higher.In contrast, for most biological indices, the correlation between W-Avg 1 and Ori tended to be higher than that between the Ori and W-Avg 1 methods.For example, Ori's BMI and the BMFI of W-Avg 1 had a 0.876 correlation coefficient, but the W-Avg 1 BMI had a 0.893 correlation coefficient, with the Ori's BMFI and the latter being relatively higher.However, this was only a relative comparison, and there was no significant difference in the correlation values; all significant differences were the same at p < 0.01.
In the fit for the model derived from the O/E derivation method, the number of species and families showed the lowest NRMSE for Ori and the highest fit; however, there was no significant difference from the NRMSE derived by the W-Avg 1 and 2 methods (Table 7).For the biological indices, the fit was higher in W-Avg 1 than in Ori, and in TESB, the difference in NRMSE between Ori and W-Avg 1 was very large.Regardless of the derivation method, the metrics with a relatively high fit were BMI, AESB, BMFI, and NRMSE, which tended to be high for the number of taxa and TESB.

Distribution of the observed, predicted, and O/E values by derivation method
The distribution of the predicted, observed, and O/E values for each metric were confirmed by applying the Ori, W-Avg 1, and W-Avg 2 methods to each reference group, and the results showed that there was a difference in the predicted and observed values between Ori, W-Avg 1, and W-Avg 2 in the number of taxa, when using the indicator taxa of BMI and BMFI, or all taxa that occurred, but no difference in the O/E values.However, for the biological indices, the differences tended to be larger between Ori and W-Avg 1 according to each derivation method, when compared with the taxon for each index.For the number of taxa, the differences between the predicted and observed values in the method using taxa rather than the derivation method meant that it could lead to a relatively more sensitive assessment, as it stemmed from the diversity of the taxa used (Cao et al. 1998;Davy-Bowker et al. 2005).However, there was little difference between the reference groups, suggesting that the proposed index taxa of the BMI and BMFI reflected most of the taxa used in the analysis, such that the sensitivity loss of the assessment was not large (Marchant et al. 1995;Smith et al. 1999).There was no difference in the O/E values, unlike the predicted and observed values, indicating that there was no difference between the Ori and W-Avg methods when the raw data were used without processing (Clarke et al. 2011).There was little difference in the distribution of the number of taxa in the lognormal distribution-based probability density function, and one-way ANOVA was used to confirm the distribution of the O/E values for each method and the differences.However, at the species level, the diversity of the taxa used for predictions and observations was higher, resulting in a higher variance at the species level than at the family level (Hawkins et al. 2000;Min and Kong 2021).
For the biological indices, in contrast to the number of taxa, there were differences in the predicted, observed, and O/E values depending on the derivation method.The deviation tended to be smaller for Ori than for W-Avg 1 in each reference group.Although the probability of occurrence for each taxon in Ori may be different, overall taxa were expected, in terms of the individual abundance, the predicted value did not show a large difference from the maximum or minimum of mean abundance for the reference sites, unlike the observed value, and the deviation was also small, and it was judged to be the result that specific taxa dominate at most sites.Furthermore, at most reference sites, the relative abundance (h i ) was also high for specific taxa such as Planariidae, Pleuroceridae, Baetidae, Isonychiidae, Heptageniidae, Leptophlebiidae, Potamanthidae, Ephemeridae, Ephemerellidae, Gomphidae, Nemouridae, Perlidae, Chloroperlidae, Elmidae, Psephenidae, Simuliidae, Chironomidae, Rhyacophilidae, Stenopsychidae, Hydropsychidae, and Apataniidae.Unlike the number of taxa, the differences in the O/E values by method were considered to be a result of the process used to drive the index, the individual abundance that was processed for the ranking, and the derivation method used for the observed values that was fixed to one method; therefore, the value for each reference group varied according to the predicted values.
For the overall reference group, among the BMI, BMSI, AESB, and BMFI, which showed higher predicted values in W-Avg 1 than in Ori, the differences between the derivation methods were not relatively large, and a statistically significant level based on the t-test was not identified in BMI and AESB.This is because the predicted abundance of taxa that increased the index value of the BMI was higher than that of the other taxa in Ori, and it was judged that the probability of occurrence for taxa with high environmental quality scores (Q i ) was high in AESB, which does not use individual abundance.Among the abovementioned families with a high relative abundance of the predicted individual abundances, most of the taxa (clear-water organisms) with a saprobic value (s i ) in the range of 0 to 1 were occupied by Isonychia, Cinygmula, Ecdyonurus, Rhithrogena, Paraleptophlebia, Cincticostella, Drunella, Isoperla, Sweltsa, Rhyacophila, Hydropsyche, and Neophylax.These taxa also had high environmental quality scores, which increased their AESB.For the W-Avg 1, as shown in the observed values, most of the reference sites were graded as A or higher; therefore, it was considered that the BMI and AESB values were high for all groups.The probability density function for O/E also showed a slight difference, only in terms of variance.
In contrast to the BMI and AESB, the BMSI and BMFI showed a relatively large difference between Ori and W-Avg 1 in their predicted values, and thus, a statistically significant difference.This was because the individual abundances of the taxa with high BMSI and BMFI values were high in the observed values, whereas the predicted abundances for taxa with lithophilic values (l i ) or medium-level sensitivity values (F i ) were high in the predicted results.Representatively, Corbicula, Potamanthus, Caenis, Macromia, and Gumaga were found to have low saprobic values for deriving BMI and they tended to be distributed in diverse substrates or fine streambeds; therefore, it was judged that the BMSI predicted by the Ori method tended to be lower than the W-Avg 1, reflecting the observed value due to the influence of these taxa.As BMFI is an index at the family level, even if there were taxa with low saprobic values at the species level, the sensitivity value was relatively low, as it was integrated at the family level (Hawkins et al. 2000;Hargett et al. 2007).Representatively, Baetidae, Caenidae, Tipulidae, Chironomidae, Hydropsychidae, and Leptoceridae belong to this category, and it was thought that there were relatively more cases with lower assessment grades by the reference group in the BMFI than in the BMI.Even in the probability density function based on the O/E value, BMSI and BMFI showed differences in their mean, median, mode, and variance between the two derivation methods; therefore, Ori was biased toward showing a larger positive value than W-Avg 1.
Unlike the above four indices, the KTI and TESB showed a statistically significant difference by showing a higher predictive value for Ori than for W-Avg 1 for all reference groups.This is also considered to be due to the high relative abundance of taxa with a relatively low thermal value (t i ; cold-water organisms), although the predicted abundance was not very low compared to the observed value.Semisulcospira, Acentrella, Baetiella, Baetis, Nigrobaetis, Ecdyonurus, Epeorus, Ephemerella, Neoperla, Simulium, and Rhyacophila were also used.Compared with KTI, the TESB showed a very large difference between the Ori and W-Avg 1 in terms of predicted value.This was attributed to the fact that although the probability of occurrence was low at all reference sites, the occurrence of most taxa was as expected, and there were many taxa with high environmental quality scores.In the probability density function based on the O/E value, the KTI showed a slight difference in the mean, median, and mode between the two derivation methods, but there was no difference in variance.In contrast, there were large differences in the TESB showed between the mean, median, mode, and variance; in particular, the variance of the O/E distribution based on W-Avg 1 was higher than the species level for the number of taxa.This was because the BMI, BMSI, and KTI were limited to a maximum of 100, BMFI to a maximum of 10, and AESB to a maximum of 5.However, the TESB was affected by the number of taxa that occurred, so there was no limit to the index value.For the taxa, a value of 1 increased steadily when an organism was added, but for TESB, a maximum of 5 was added, and the change was larger than the number of taxa.Therefore, it was considered to have the highest degree of variance when compared with the other indices (Min and Kong 2021).
Unlike the number of taxa, where there was no difference between the derivation method and the method for using the taxon in the O/E of each reference group, the biological indices showed a difference between the O/E and the derivation method according to the differences in the observed and predicted values, and the tendency was different for each index; however, for most groups, the W-Avg 1 method resulted in the O/E being close to 1.0.This is because the W-Avg 1 derives a predicted value based on the observed results.The Ori also expects the predicted abundance based on actually observed individual abundances, and derives biological indices based on this, however at this time, as mentioned above, there was no significant difference in the probability of occurrence of taxa by reference sites, or the ranking of relative abundance based on predicted abundance.For this reason, for the Ori, there were cases where a low index value was derived or an excessively high index value was obtained in a relatively good environment, and this was very prominent in the TESB.

Degree of overlap and fit between indices based on Ori and W-Avg 1
For the reference sites, the O/E value for each derivation method showed a Pearson's correlation coefficient of ≥ 0.8 between the same metrics and was statistically significant at < 0.01.At this time, the number of taxa, TESB, BMI, AESB, and BMFI showed high levels of correlation as they were indices with similar characteristics.KTI, which represents water temperature according to altitude, and BMSI, which represents the type of streambed based on substrate composition, showed relatively high correlations with BMI, AESB, and BMFI, which represent organic pollution, compared with other indices.
If the habitat or microhabitat of organisms in the river becomes simpler owing to disturbance, the diversity of taxa also decreases, which can lead to the deterioration of the river ecosystem (USEPA 1999).A biological community is a receptor that reflects the current state of various environmental variables in a complex manner without being biased toward either the physical or chemical variables of the river (Urrea-Clos and Sabater 2012).Environmental variables that determine the distribution of biological communities can interfere with each other.For example, changes in streambeds due to sedimentation caused by a decrease in discharge and velocity can cause the deterioration of water quality (Judy et al. 1984;Klemm et al. 1993).Conversely, an increase in organic and inorganic pollution causes degradation of the streambed, and an increase in nutrients causes the flourishing of attached algae on the substrate (Mebane et al. 2021).It is thus important to analyze the distributions of biological communities by utilizing the complexity of chemical and physical variables; however, it is important to understand environmental changes (Yoon et al. 2018), as it is difficult to distinguish the effects of physical and chemical variables on the distribution of benthic macroinvertebrate communities.In view of the correlation between the indices, the water quality and physical variables in this study acted in a complex manner rather than having a distinct effect on the distribution of biological communities.However, considering only the reference sites with the premise that most environmental variables were good, rather than the data of all sites used in this study, this was judged to be a result entirely derived from the river continuum concept proposed by Vannote et al. (1980).In other words, when a river flows from upstream to downstream, the environmental variables gradually change.In general, when the stream size increases, the water temperature decreases alongside the altitude, the velocity at the sites (or points) accessible to people slows, water depths deepen, streambeds become fine due to sedimentation, and water pollution tends to increase.Considering these characteristics of rivers, taxa living upstream prefer clean water quality, coarse streambeds, and low water temperature, whereas downstream taxa are relatively resistant to pollution and have adapted to various river substrates and water temperatures.In the middle stream, the taxa between the upper and lower streams were redundant.This suggests that even if an index represents different environments, one index can indirectly represent different environmental variables.
In terms of the fit of the model for each derivation method, the difference in NRMSE was insignificant, as it was< 0.02 for the number of taxa, BMi, and AESB, and the difference in BMSI, KTI, TESB, and BMFI was approximately 0.1, showing relative but not large numerical differences compared to the former case.In contrast, TESB showed a very large difference of 7.514.The fact that W-Avg 1 exhibited the largest NRMSE for the number of taxa was attributed to the range of taxa used.Conversely, for all biological indices, W-Avg 1 showed a lower NRMSE than Ori; thus, the predicted values were closer to the observed values.This is because W-Avg 1 is derived from the the observed value of the index.
In BMI, BMSI, KTI, AESB,and BMFI, where the differences in NRMSE were smallerderiving the predicted index value based on the observed index showed a higher degree of fit than the index derived based on the more predicted abundance.This shows that W-Avg 1 is more suitable than the existing method for the assessment of biological aquatic environments in the Republic of Korea.As TESB is an index that is the sum of the environmental quality scores for each taxon that has occurred, all taxa with a probability of occurrence were used.Therefore, the TESB of Ori was the same for each reference site and was predicted to be very high, making it impossible to compare with the observed value.This is also the best indication that the newly proposed prediction method is more suitable for the assessment of rivers in the Republic of Korea.
The purpose of this study was to determine a new predicted-value derivation method for the Korean River Macroinvertebrate Prediction and Assessment System (KRIMAPS), which was proposed to assess the integrity of rivers in the Republic of Korea based on long-term cumulative data, and to generate a more accurate and high-precision assessment method by comparing existing methods.There was no difference in the number of taxa when the distributions of the predicted, observed, and O/E values derived from Ori and W-Avg 1 or 2 were compared, and the correlation and fit of the O/E for each derivation method analyzed.However, for biological indices, the new proposed method (W-Avg) for deriving the predicted index value based on the observed index value was more suitable for assessing rivers in the Republic of Korea than the existing method (Ori) for predicting the index value based on the predicted abundance of each taxon.Therefore, this study proposes W-Avg 2 as it has less deviation, and only the indicator taxa of BMI or BMFI are required for national biological monitoring (such as NAEMP).As all taxa were used for the detailed surveys conducted on the specific streams, W-Avg 1 was also found to be favorable as it has a relatively high assessment sensitivity.
In this study, the number of taxa and biological indices based on benthic macroinvertebrates developed in the Republic of Korea were analyzed as the main targets; however, additional examination of community structure metrics, such as dominance (McNaughton 1967), diversity (Shannon and Weaver 1949), richness (Margalef 1958), and evenness (Pielou 1975) will be required in the future.It will also be necessary to review the applicability of the method and to propose a assessment grade.

Figure 1 .
Figure 1.locations of the 357 sampling sites (gray circles) in the republic of Korea.thick black lines indicate the watersheds of the four major rivers (han, nakdong, geum, and yeongsan-seomjin).thin blue lines indicate the main, tributary, and other independent streams within each watershed.

Figure 2 .
Figure 2. two-way indicator species analysis (twInsPan) of the 357 reference sites based on the benthic macroinvertebrate communities.the level represents twInsPan divisions and the numbers in parentheses indicate the number of reference sites for each group (Min and Kong 2021).

Figure 3 .
Figure 3. Distribution patterns of each environmental variable based on the two-way indicator species analysis (twInsPan) groups for the 357 reference sites, which were assessed using the benthic macroinvertebrate communities (a) catchment area, (b) altitude, (c) latitude, (d) longitude, (e) water depth, (f) velocity, (g) mean diameter of streambed, (h) 5-day biochemical oxygen demand (BoD 5 ), and (i) total phosphorus (tP).the lines on the bar graph indicate standard deviation.

Table 2 .N
Metrics used by the two categories for the derivation of predicted (p), observed (o), and observed to expected (o/e) values.Biological index Benthic Macroinvertebrates Index (BMI;

Figure 4 .
Figure 4. Distribution of the observed, predicted, and ratio of observed to predicted (o/e) values for the number of taxa at the species-level (a-c) and family-level (d-f), Benthic macroinvertebrates index (BMI, g-i), and Benthic macroinvertebrates streambed index (BMsI, j-i) for each reference (twInsPan) group.the original method (ori), weighted average method 1 (w-avg 1), and weighted average method 2 (w-avg 2) refer to the taxa and derivation methods used to calculate observed, predicted, and o/e values, respectively.each letter and the corresponding color-coded line indicate the evaluation class interval for each biological index.Below the D class line indicates e class.

Figure 5 .
Figure 5. Distribution of the observed, predicted, and ratio of observed to predicted (o/e) values for the Korean thermality index (KtI, a-c), total ecological score of Benthic macroinvertebrate community (tesB, d-f), average ecological score of Benthic macroinvertebrate community (aesB, g-i), and Benthic macroinvertebrates family-level biotic index (BMfI, j-i) for each reference (twInsPan) group.the original method (ori) and weighted average method 1 (w-avg 1) refer to the taxa and derivation methods used to calculate the observed, predicted, and o/e values, respectively.each letter and the corresponding color-coded line indicate the evaluation class interval for each biological index.Below the D class line indicates e class.

Figure 6 .
Figure 6.Distribution of the observed to predicted (o/e) values for the number of taxa as the species-level (a) and family-level (b), Benthic macroinvertebrates index (BMI, c), Benthic macroinvertebrates streambed index (BMsI, d), Korean thermality index (KtI, e), average ecological score of Benthic macroinvertebrate community (aesB, f), total ecological score of Benthic macroinvertebrate community (tesB, g), and Benthic macroinvertebrates family-level biotic index (BMfI, h) for the 357 reference sites.the original method (ori), weighted average method 1 (w-avg 1), and weighted average method 2 (w-avg 2) refer to the taxa and derivation methods used to calculate o/e values, respectively.

Table 1 .
reference conditions of each environmental variables by stream types considered for selecting of reference sites.

Table 3 .
Mean, median, mode, and standard deviation of the predicted to observed (o/e) values based on the probability density function for the lognormal distribution for each derivation method.

Table 5 .
T-test results for expected to observed (o/e) of each biological index.Korean thermality Index, tesB: total ecological score of Benthic macroinvertebrate community, aesB: average ecological score of Benthic macroinvertebrate community, BMfI: Benthic Macroinvertebrates family-level Biotic Index when the p-value of leven's equality of variance test was higher than 0.05 and the pooled variance estimate lower than 0.05, the Behrens-fisher statistic was applied for the t-test p-value.

Table 4 .
one-way anoVa results for the expected to observed (o/e) of number of taxa at the species and family levels.

Table 7 .
normalized root mean square error (nrMse) of the observed to predicted (o/e) values derived from original method (ori), weighted average method 1 (w-avg 1), and weighted average method 2 (w-avg 2) for each metric.

Environmental quality score (Kong, Park, et al. 2018) for each taxon to derive total ecological score of benthic macroinvertebrate community (TESB) and average ecological score of benthic macroinvertebrate community (AESB)
the higher the environmental quality score, the higher the sensitivity to organic pollution.Kub is unnamed (or undetermined) taxa.

Appendix 5. Sensitivity value for each taxon to derive benthic macroinvertebrates Family-Level biotic index (BMFI; Kong et al. 2019)
(Continued)