Assessment of medium and small river health based on macroinvertebrates habitat suitability curves: a case study in a tributary of Yangtze River, China Free

River ecosystems are an important component in natural ecosystems. They link terrestrial to aquatic ecosystems in terms of exchanging substances, energy, and information (Xiao et al., 2019; Yu et al., 2019). River ecosystems have been disturbed by hydropower development and human activities (e.g., urbanization), resulting in deteriorating water quality and loss of riverine biodiversity (Silva et al., 2014; Moosová et al., 2019). At a global scale, this seriously hinders the sustainable development of society and economy. It is thus important to evaluate the health of river ecology before effective restoration practices can be made (Boulton, 1999).

After proposing the concept of ‘river health’, Europe initially conducted river health assessment that focused on physical and chemical indicators (Patrick, 1973). Later the composition and structure of aquatic organisms including phytoplankton, zooplankton, benthic invertebrates and fish were introduced into the evaluation of river health (Smith et al., 1999; Zhang et al., 2018). During the past decades, a large number of assessment methods have been developed and they can be in general assigned into two groups: (1) comprehensive indicator method and (2) biological monitoring method. The comprehensive indicator method was first developed in Europe in the 1900s with a primary focus on flow velocity, water depth and water quality (Patrick, 1973), e.g. Index of Stream Condition (ISC) (Ladson et al., 1999). It has evolved to the current combined assessment method in which chemical, physical, biological parameters, and human activities in the river ecosystems are all considered (Luo et al., 2018; Lin et al., 2019). The biological monitoring method attempts to quantitatively describe biological integrity by integrating biome composition, structure, species traits and functional parameters. The earliest applications of the indicator species method for river health assessment were River Invertebrate Prediction and Classification System (RIVPACS) (Wright et al., 1984) and Australian River Assessment System (AUSRIVAS) (Sudaryanti et al., 2001). In recent years, methods such as Multi-Metric Indices (MMI), have become increasingly popular, and are widely used to evaluate the health of large rivers (Chen et al., 2014).

The comprehensive indicator method and biological monitoring method achieved a success in major rivers. In medium and small rivers, however, they cannot be readily applied. Their distinct characteristics like short channel length, small watershed area, and distributed in remote areas make it difficult to collect enough environmental variables and other biological descriptors to establish a comprehensive evaluation system. The health status of medium and small rivers has become a global concern. The EU Water Framework Directive (WFD) (Kallis & Butler, 2001) has been implemented in all EU member states since 2000 and has been developed to investigate study methods and assessments of water status for virtually all European watersheds. The USA and Canada have used predictive models, e.g., Watershed Health Assessment Tools-Investigating Fisheries (WHAT-IF) and Pressure-State-Response (PSR) (Bricker et al., 2003; Rashleigh et al., 2006), which consist of hydrology, chemistry and biology variables to diagnose the health of Muskoka and Mississippi watershed.

Assessing the health of medium and small rivers is important not only for China, but also for all other countries. As a result of lack in sufficient monitoring data for Chinese medium and small rivers, the applications of evaluation methods that were developed by Europe and USA/Canada have been limited in China. To fill this gap, the improved method, macroinvertebrates habitat suitability curves (M-HSC) was proposed in this study (Figure 1). Macroinvertebrates are an important component of benthic food chain in the medium and small rivers. They play a vital role in material cycling and energy flow of river ecosystems. Their community composition and distribution are sensitive to aquatic environmental changes, which makes it a suitable indicating organism for river health assessment. The macroinvertebrates density can be an indicator of the health status of rivers.

The habitat suitability index of aquatic organisms is a quantitative evaluation method that describes the relationship between habitat preference and habitat factors of aquatic organisms. The core of this method is to build a continuous suitability curve of the relationship between the preferences of target species for individual habitat factor. The suitability curves use a value between 0 and 1 to define the preference of the target species for a single factor. The peak of the curve represents the optimal or favorite range for the factor takes the value of the habitat factor as the abscissa and the preference degree of the target species to the habitat factor (i.e., suitability, SI) as the ordinate (Armour & Taylor, 1991). Common methods for constructing suitability curves include generalized additive models (GAMs), generalized linear models (GLMs), fuzzy logic models, artificial neural networks (ANNs) and classification trees. To compare the performance of these models, researchers (Gevrey et al., 2005) determined functional feeding groups and evaluated species richness for Dutch macroinvertebrates using five modeling techniques (ANN, Regression trees, GAM, Partial Least Squares, and Multiple Linear Regression), which showed that the ANNs and GAMs performed better concerning to root mean square errors and correlation coefficients compared with the other three methods. Due to its ability to deal with highly nonlinear, non-monotonic response curves, GAMs have been used to predict the correlation between species distribution and the environment such as the distribution of fish, phytoplankton, and benthic macroinvertebrates in the oceans, estuaries, rivers and lakes (Mouton et al., 2011; Yi et al., 2014).

The purposes of this study are: (1) to enhance the accuracy of macroinvertebrates habitat suitability curves based on Generalized Additive Model, (2) to propose a new river health assessment method to assess medium and small river health based on macroinvertebrates habitat suitability curves (M-HSC), and (3) to apply this new method to evaluate the health condition of Qiaobian River, a tributary of Yangtze River, China. The successful development of this assessment method provides new insights for future river health assessments.

Qiaobian River (QBR), a first-order tributary on the right bank, is located in the mid-lower reach of the Yangtze River, China (Figure 2(a) and 2(b)). QBR has a watershed area of 295 km² and a total length of 40 km, is one of the important drinking water sources in Yichang. The highest elevation in QBR basin is 1,168 m, with an average hydraulic slope of 11.4‰. Three villages are distributed in the QBR basin. They are Dianjun, Qiaobian and Tucheng. The total population in the QBR basin is 39,969. QBR is an important source of drinking water in Yichang, which plays an essential role in protection against floods and in providing habitat for aquatic organisms. Prior to the 2010s, the QBR basin is one of the most famous tourist attractions in Yichang city due to its well conserved nature. Since the 2010s there has been an increase in water pollution due to urbanization, industrial development and climate change, which poses a serious threat to drinking water safety. In this study, a total of 8 sampling sites representing a range of different habitats were sampled for benthic macroinvertebrates and ecological variables in winter 2018 and summer 2019 (high and low flow seasons, respectively) (Figure 2(c)).

Three samples were collected from the wadable shallow area within 100 m stretch of each sampling site using a 500 μm D-shape net. Six subsamples (0.3 m² each) were collected at each sampling site with a total area of 1.8 m² investigated. The macroinvertebrate samples were sieved through a steel mesh (pore size 0.5 mm) and then transferred into 500 mL bottles. The samples were fixed with 10% formaldehyde solution for storage. Manual picking and identification were carried out indoor on a white ceramic tray (5 cm × 5 cm grid size). Identification of all individuals was done to the lowest possible taxonomic unit, in which oligochaetes were identified to the class and chironomids to subfamilies. Density of macroinvertebrates was calculated as the number of individuals per area (ind m⁻²). Only the head was counted when the samples were damaged. Biomass of macroinvertebrates (g m⁻²) was obtained by weighing the samples. Filter papers were used to dry the surface of the samples before weighing.

Seasonal field sampling was conducted in December 2018 and April 2019 by a group of 5 people. All the field sites were sampled within 3 days in each season. Water depth (Dep), water temperature (Temp), turbidity (Turb), pH and dissolved oxygen (DO) were measured in situ using a multiparameter probe (Hydrolab DS5, Hach, USA). Flow velocity (Vel) was measured using an Acoustic Doppler Velocimetry (ADV, Nortek, Norway). Median particle size (D₅₀) of riverbed sediment was measured using a laser diffraction particle analyzer (Microtrac S3500, Microtrac, USA). Water samples were collected at each sampling site kept in a cold box in darkness and transported back to laboratory for water chemistry analysis. Parameters including chemical oxygen demand (COD_Mn), total nitrogen (TN), total phosphorus (TP), ammonium (NH₄⁺), nitrate-nitrogen (NO₃⁻) and phosphate (PO₄³⁻) were determined.

To avoid subjective factors, the weight of each suitability index is calculated using entropy method. There are mainly three steps to calculate the weight with entropy method.

• 1.

  Normalization of raw data matrix.

• 2.

  Define entropy

• 3.

  Define entropy weight

A total of 1,152 benthic macroinvertebrates were collected, which could be assigned to 3 departments, 5 classes and 19 families (Table 1). Among them, arthropods were the most abundant species, accounting for 68.4% of the total, followed by mollusks, accounting for 26.3% of the total species, and the least number of animals, accounting for 5.3% of the total species. The main dominant species were Chironomus flaviplumus, Ilyoplax deschampsi, Semisulcospira cancellate, Unio douglasiae, Radix auricularia, Bellamya purificata and Corbicula fluminea (Table 2). Corbicula fluminea was chosen as an indicator species here because it was found in all sampling sites.

QBR is a typical shallow river (Table 3, max. water depth 3.1 m). In general, QBR was in a low level of pollution that was characterized by low nutrient and COD_MN. The spatial variabilities of these variables were not significant except for sediment size. The upstream was dominated by gravel (D₅₀ = ∼50 mm), in contrast to silty downstream (D₅₀ = 0.014 mm).

A CCA was done to explore correlations between benthic macroinvertebrate data and chemical variables (Figure 3). The first two axes explained 14.75 and 13.22% of the variation, respectively. A total of 47.72% of the variation was explained by the first four axes (p < 0.05). These water chemistry parameters including TN, NO₃⁻, NH₄⁺, COD_Mn, pH, and DO had significant effects on axis 2 (Figure 3, Table 4), in which TN, NO₃⁻, NH₄⁺, COD_Mn were the variables with the greatest negative impact and pH, DO had the greatest positive influence. The angle between the dominant species, Corbicula fluminea, and DO was acute, indicating that Corbicula fluminea positively responded to DO; the angle between Corbicula fluminea and TN, NO₃⁻, NH₄⁺, COD_Mn and pH was obtuse, suggesting a negative response of Corbicula fluminea to these parameters.

A CCA was done to explore correlations between the benthic macroinvertebrate data and chemical variables (Figure 4). The first two axes explained 23.28 and 20.04% of the variation, respectively and a total of 67.40% was explained by the first four axes (p < 0.05). D₅₀, water depth, and Turb had a significant negative influence on axis 2 (Figure 4, Table 5). The obtuse angle between the dominant species, Corbicula fluminea, and D₅₀, water depth, and Turb suggested a negative response of Corbicula fluminea to these physical parameters.

Pearson correlation coefficient was done to eliminate highly correlated variables and 0.75 was used as the threshold (Table 6). Among the chemical factors, a high correlation was found between COD_Mn and NH₄⁺, TN and NO₃⁻, TN and NH₄⁺, respectively. Based on data availability, NH₄⁺ and NO₃⁻ are excluded. D₅₀ and Temp were also excluded due to the higher correlations with other variables. Then the main environmental parameters driving the distribution of macroinvertebrates were COD_Mn, TN, pH, DO, Turb and water depth.

A GAM was set up to quantify correlations between the habitat suitability of Corbicula fluminea and key environmental variables including COD_Mn, TN, pH, DO, Turb and water depth (Figure 5).

According to the definition and calculation method of entropy weight method, the calculated results of habitat suitability factor weight of macroinvertebrate in QBR are shown in Table 7. The weights of each indicator were sorted as Dep > TN > pH > DO > COD_Mn > Turb.

There are six habitat indicators and the highest value of each is 1. The total score was divided into five equal parts, and the river health evaluation standards were constructed according to the classification of Nature, Health, Sub-health, Unhealth and Ill (Table 8).

Through the interpolation calculation, the distribution of the health index in QBR was obtained (Figure 6). In general, the health status of QBR declined streamwise. The upstream (site S7-8) was in a good health and could be assigned into Class I; this is followed by a gentle decrease at site S6 which fell in Class II; a sharp decrease in evaluation score resulted in downgraded middle reaches to a sub-health status with site S5 in Class III and site S3-4 in Class IV; the most downstream sites S1-2 fell into class V, suggesting a bad health status.

The CCA ranking results reveal the composition of habitat factors that has an important impact on the health of macroinvertebrates' habitat in QBR. The quality of macroinvertebrates' habitat is determined by water chemistry (TN, NO₃⁻, NH₄⁺, COD_Mn and pH) and physical parameters (D₅₀, Dep and Turb). This is consistent with a large number of previous studies, i.e., in riverine ecosystems the health of macroinvertebrates' habitats is affected by a variety of physical and chemical factors (Karr & Chu, 2000; Freeman & Rogers, 2003). Pearson correlation analysis can eliminate the more relevant factors and reduce the redundancy of indicator library. Based on the analysis, a total of six indicators (COD_Mn, TN, pH, DO, Turd and Dep) were selected as the key habitat factors affecting the habitat health of dominant species of macroinvertebrates in QBR. The results were similar to other studies (Cai et al., 2017; Yang et al., 2018). While factors such as NH₄⁺ and D₅₀ were excluded in this study, these factors still play an important role in the health of macroinvertebrates' habitat (Hilsenhoff, 1988; Yi et al., 2016a).

The most dominant species in the entire river is Corbicula fluminea. Besides the dominance of Corbicula fluminea in the upper reaches, Ilyoplax deschampsi was another dominant species; in the downstream, Chironomus flaviplumus were also one of the major dominant species. The succession of dominant species of macroinvertebrates is closely related to changes in river water quality. The water quality in the upper reaches of QBR was relatively good, making it an ideal habitat for less pollution-resistant species such as Ilyoplax deschampsi. The lower reaches of QBR flow through densely populated areas (Qiaobian town and Wulong community), which led to deteriorated riverine water quality. This was clearly indicated by the apparent dominance of Chironomus flaviplumus, a pollution-resistant species.

A reasonable habitat suitability model is of great importance to river restoration and management (Yi et al., 2016b). Currently there are many methods for estimating habitat suitability of aquatic organisms such as preference curve method (Bockelmann et al., 2004) based on experts' experience, generalized linear model and GAM (Vasconcelos et al., 2013; Yi et al., 2017) etc. The comparison of performances of different models showed that GAM model was an ideal method for habitat suitability analysis because of its flexibility, ability to simulate complex nonlinear relationships, and simple and accurate calculation (Yi et al., 2018).

This study used GAM to explore the mutual selection relationship between dominant species of macroinvertebrates and key habitat factors in QBR – habitat suitability. The results suggested that the habitat suitability of the Corbicula fluminea in QBR was negatively correlated with COD_Mn. The habitat suitability decreased to 0 when COD_Mn concentration exceeded 3.5 mg L⁻¹. The optimal TN concentration of Corbicula fluminea in QBR is 0.27 mg L⁻¹, and the increase of TN concentration leads to a sharp decline of habitat suitability. As Corbicula fluminea is a species with moderate resistance excessive nutrients make them difficult to survive (Burket et al., 2019; Siuda et al., 2020). In a normal range, the habitat suitability of Corbicula fluminea was well correlated to DO with a positive response (Riding et al., 2019), but supersaturated DO could also inhibit the spontaneous breathing process of aquatic organisms, leading to its death (Shen et al., 2019). This study showed that when DO was >13 mg L⁻¹, the habitat suitability of Corbicula fluminea was reduced to 0. The impact of DO oversaturation on macrobenthos has rarely been reported, which is worth investigating in the future. The optimal water depth for Corbicula fluminea in QBR was 0.35 m, and the habitat quality was highest when the water depth was 0.16–0.5 m (Beauger et al., 2006). The most suitable pH for Corbicula fluminea is 8.41, which was consistent with the optimal pH for fish (Leivestad & Muniz, 1976). The most suitable turbidity in QBR was 1.76 NTU, and Turb was negatively correlated to the habitat suitability of Corbicula fluminea. However, the habitat suitability of Corbicula fluminea appeared to rise briefly when Turb was in the range of 4–6 NTU. This might be related to the release of nutrients associated with disturbance to the riverbed to provide food for Corbicula fluminea (Müller et al., 2019).

The Habitat suitability value of different habitat indexes for macroinvertebrates and fishes (common carp, Spinibarbus hollandi and so on) (Edwards & Twomey, 1982; Li et al., 2015) showed similar results (Table 9). The common carp (Cyprinus carpio) is not only common in Asia, it is now also found on many continents, whose habitat status becomes a global concern. Fishes' habitat suitability value (12 mg L⁻¹, 0.20–0.60 m, 6.00–10.00, 10.60 NTU), the value of DO, Dep, pH, Turb, is the same with macroinvertebrates' (11.17 mg L⁻¹, 0.35 m, 8.42, 10.16 NTU). As for the value of COD_Mn and TN, Fishes' is a little higher than macroinvertebrates', most likely because large aquatic species like fishes can absorb and utilize more nutrients in water column than smaller organisms. However, low nutrients content indicates less environmental pollution and healthier habitat conditions, which is suitable for species at high-trophic states. The results based on macroinvertebrates could be extrapolated to predict the large species of interest such as fish.

The water quality parameters of QBR in each sampling season are shown in Figure 6. Little changes in pH were found between each site and sampling seasons. In contrast, COD_Mn, DO, and TN showed different trends. From S8–S1 (upstream to the estuary), the concentration of COD_Mn and TN showed an increasing trend, which is different from DO that decreased streamwise. These three factors all responded negatively to the increase of flow, indicating a clear dilution effect.

The river was spatially heterogeneous in river health status (Figure 6). The health status of upstream was significantly better than that of downstream with a clear declining trend from upstream to downstream. This can be explained by the higher flow velocity and fewer inhabitants at the upper reaches of QBR than at the downstream. The sudden decrease of health status in the middle reaches might be partly due to the drainage of sewage water from the surrounding towns (Tucheng and Qiaobian). The heavily deteriorated health at lower reaches of QBR can be attributed to the development of industry in this area, where untreated industrial wastewater was discharged to the river.

The water quality (represented by the concentration of COD_Mn, DO and TN) was also temporally variable in QBR (Figure 7). These environment variables were higher in low flow season than in high flow season. This can be explained by the higher anthropogenic pollution at downstream reaches of QBR than at the upper reaches. The population in the downstream area, the Dianjun community, is nearly three times than in Qiaobian and Tucheng town. All domestic sewage waters from cities and towns are discharged into rivers through sewage pipes. The excessive urban population exceeds the treatment capacity of the local wastewater treatment plants and the sewage discharged from the sewage outlet of the Dianjun contains high nutrients. In the low flow season, as the smaller river flow and velocity, the sewage discharged from the sewage outlet cannot be diluted and flushed away by the river runoff, resulting in the poor water quality of QBR in the low flow season. It is suggested that the Yichang Municipal Government should take an acute measure during the dry season to reduce the pressure brought by human activities in QBR.

River health assessment for QBR was performed by Yichang Environmental Protection Agency (YEPA) using the environmental quality standards for surface water in China, GB3838-2002 (SEPA & AQSIQ, 2002). The results showed a good status of habitats in QBR, and the average values of water chemistry such as TN and TP (0.49 and 0.05 mg L⁻¹, respectively) did not exceed their specified limits (0.50 and 0.10 mg L⁻¹). Only the habitats at lower reaches of QBR were in poor condition with TN (1.06 mg L⁻¹), NO₃⁻ (0.70 mg L⁻¹) and TP (0.17 mg L⁻¹) exceeded the prescribed thresholds. Our assessment results were in a good agreement with the data from YEPA (Figure 8). This suggests that assessment based on habitat suitability of macroinvertebrates can well reflect the actual situation of river health.

While our study is in line with field monitoring results, certain limitations remain for the health assessment by considering population (dominate species) only. Future work should include the correlation between species abundance and river health. Moreover, it is noteworthy that the present study did not consider the complexities of food web. Food chain parameters such as aquatic algae and organic matter content/quality should be included in the future habitat suitability model.

This study proposes a new river health assessment method based on macroinvertebrates habitat suitability curves, standard curves method (SCM), to assess health status of medium and small rivers. In this paper, as a case study, we investigated the distribution of macroinvertebrates, chemical and physical parameters during winter 2018 (low flow season) and summer 2019 (high flow season) of Qiaobian river, a tributary of Yangtze River, China. Six variables including COD_Mn, TN, pH, DO, Turbidity and water depth were used to describe the river health suitability of macroinvertebrates. Meanwhile, Generalized Additive Model (GAM) was applied to improve the accuracy of suitability curves. The results showed the health status of QBR is spatially heterogeneous being apparently better upstream than downstream.

Compared with the environmental quality standards for surface water in China (GB3838-2002), the new river health assessment method based on macroinvertebrates habitat suitability curves (M-HSC) showed good uniformity and consistency. This study is a new attempt to assess medium and small river health and provide insights for future river health assessments.

This study was financially supported by the National Natural Science Foundation of China (grant No. 51979149). We thank Yichang Environmental Protection Agency (YEPA) for providing data.

The authors declare no conflicts of interest.