Assessing limitation factors and thresholds for macroinvertebrate communities in response to land use gradients

. Transformations of land use from natural to anthropic type have been recognized as a significant trigger which degenerate the aquatic ecological quality seriously. However, there was still lack of enough evidence which the extent of changes in land use should be set as a biodiversity conservation target to protect aquatic ecosystem. To understand the corresponding variations of aquatic organisms to environmental gradients and set the conservation threshold values for land use, data of physicochemical parameters and macroinvertebrate communities were sampled in the Hun-Tai River Basin during 2009 and 2010. The main objectives of the present study were (i) to explore limiting factors that affect the distribution of macroinvertebrate communities with land use gradients, (ii) to estimate thresholds for the conservation of macroinvertebrate communities derived from generalized additive models (GAMs) and Threshold Indicator Taxa ANalysis (TITAN), respectively. The results indicated that macroinvertebrate communities’ structure and integrity were strongly negative with nutrient, organic contaminants content, %CropArea and %ImperviousArea. Under a precaution perspective and given current levels of land use, this research might provide some useful strategies for appropriate land exploitation management and improving water quality and biodiversity conservation in river ecological restoration.


Introduction
Changes in land use were regarded as the most severe threats to global biodiversity conservation [1,2] . Some associated affection were modifications of biological communities [3] , and the loss of ecological processes and ecosystem services associated with food supply, freshwater and forest resources [4] . Freshwater ecosystems have become among the most threatened compared to the terrestrial and marine ecosystems [5] . This mainly due to settlements of humans near waterways were disproportionate and riparian zones were modified extensively [6] . These modifications of habitat were the most significant threats to freshwater ecosystems through some direct and indirect pathways [7] , and every 10% of lost natural catchment land use led to an average loss of 6% of native freshwater fish and macroinvertebrate species.
Currently, humans have extensively changed the Earth's land surface, around 40% of the Earth's ice-free land surface have already used for anthropogenic new habitats [8] . Several studies demonstrated that ecosystem management and/or restoration strategies have no positively significant effect on aquatic biota without consideration for land use changes in catchment scale [9,10] . Thus, it has become a major obstacle when facing protecting freshwater biodiversity.
Ecological thresholds are considered as transition points or zones of relatively rapid change that response to small, continuous changes of alternate ecosystem states or ecological condition. These abrupt transition points are particularly associated with anthropogenic environment gradients that represent novel physicochemical conditions [11] . Species distributions, under the critical novel environmental gradient, may be resulted in simultaneous and abrupt change. Threshold exploration depends on different response variable [12] , assumed shape of the response and appropriate statistical model [13] , which are closely relevant to the location of a threshold or whether a threshold exists [14] . Combining different analysis approaches to separate responses would also help managers to identify dominant negative and positive responses within assemblages, and estimate where major shifts in assemblage structure occur. However, current studies were lack of enough evidence that utilize these approaches in assessing responses in assemblage level to habitat loss which driven by changed land use.
In this study, our purpose is to identify ecological thresholds of land use for macroinvertebrate community combined different approaches. Macroinvertebrate play a critical role in the energy pathway and nutrient cycling. Despite some bioindicators, diatoms, macrophytes and fish have been used to assess the water quality [15,16] , macroinvertebrate have been the most widely used as indicators for biological assessment [17] . Thus, it is ICESD 2021 imperative to understand the critical change points for macroinvertebrate community and individual species stressed by land use and provide useful information to managers. The main purposes of this paper are as follows: (1) to explore the effects of environmental parameters on macroinvertebrate community with disturbance gradient; (2) to estimate the main water quality parameters thresholds and land use thresholds for macroinvertebrate community.

Study area
The Hunhe and Taizi River Basin (hereafter, Hun-Tai River Basin) (40°40′-42°10′ N, 122°05′-125°17′ E) is located in the southeast region of Liaoning Province, which is two of the largest tributaries of the Liaohe River. The total drainage area and the annual average discharge is 2.5×10 4 km 2 and 6.9×10 9 m 3 , respectively. The Hun-Tai River Basin is an important agricultural and residential region in Liaoning Province. The CropArea% and impervious% is 30.1% and 8.6% of the total drainage area in the Taizi River Basin, and 33.4% and 15.5% for Hunhe River Basin. During 1990~2010, the areas of agriculture and urban cover in this region indicated a significant amplification [18] . Landscape fragmentation, reduction of wetland area and increase of impermeable area, led to a degradation of ecosystem in this basin [19] .

Acquisition and identification of macroinvertebrates
To avoiding others stresses, such as point pollution driven by industry, flow dynamics variation caused by dam construction etc., 256 samples were collected away from the urbanization of commercial and industrial areas from 2009 to 2010. At wadeable sampling sites, macroinvertebrates were sampled within the 100 m reach by using Surber net (30×30 cm, 1 mm mesh size), three replicates were obtained from two riffles and one shallow pool (<40 cm deep). For non-wadeable sites, three replicates along the riverside were sampled by using Surber net at the shallow water, a D-frame net is used for the additional replicates in the deep-water areas (>40 cm deep) [20] . All the samples were preserved within 95% ethanol. In the laboratory, all the specimens are sorted, counted, and identified to the lowest taxonomic level.

Water sampling and analysis methods
Eleven parameters were measured, including water temperature (WT, o C), pH, electrical conductance (EC, μS/cm), total dissolved solids (TDS, mg/L), dissolved oxygen (DO, mg/L), suspended solids (SS, mg/L), total nitrogen (TN, mg/L), total phosphorus (TP, mg/L), 5-day biological oxygen demand (BOD 5 , mg/L), permanganate index (COD Mn , mg/L), ammonia nitrogen (NH 3 -N, mg/L). Generally, the values of WT, pH, DO, EC, and TDS were directly determined in situ using a portable multiparameter water quality monitoring instrument. Water samples were preserved in polyethylene plastic bottles and kept below 4 o C, for used the others parameters which were measured in laboratory according to methods associated with environmental quality standards for surface water [21] .

Extraction of land use
Combined with ground-truthed during field surveys, land use patterns were extracted through buffer tools in ArcGIS in autumn 2010 before the crops harvest. For avoiding the influence of great temporal variability of individual plots and limiting by solar illumination and cloud cover, multi-temporal data sets were applied to increase the classification accuracy. Besides, the Landsat images (30m resolution) were also orthorectified to avoid misclassification by using a Digital elevation model and field surveys. Seven patterns of land uses were classified: vegetated land, crop land, impervious layer, waters, unused land, wash land and wetland. The upstream with 10 km and riparian width with 1 km for each sampling site were focused and measured the land use composition.
To understanding the effect of land uses on river ecosystem, the current study focused on the three dominant land use patterns, including forest, crop and impervious layer.

Data analysis
Summary statistics were computed for the basic water quality parameters. Detrended Correspondence Analysis (DCA) was performed to choose the appropriate type of model for ordination analysis [22] . In this study, the canonical correspondence analysis (CCA) would best fit the data and explore the relationship between environmental parameters and structure of macroinvertebrate community for the value of the longest gradient length is > 4. All data were log (x+1) transformed. Species occurred no less than three or more sites and the relative abundance no less than 1% were removed from the analysis to reduce the influence of rare species [22] . Monte Carlo test with 999 permutations and the forward selection option were used to assess the contribution and significance of each environmental parameter that explained the majority of variance in the species data.
Principle Component Analysis (PCA) was performed to explore the main gradients with all the environmental parameters by using SPSS 22.0. The KMO and Bartlett's sphericity test and Eigenvalues greater than 1 were taken as criterion for the extraction of the principal components. The parameter with absolute loading value > 0.60 was considered as a strong positive one. The proportional variables were arcsine squared root transformed and the continuous variables were log(x+1) transformed to obtain normal distributions.
Generalized additive regression models (GAMs) were used to fit continuous response relationships of macroinvertebrate community integrity to each stress gradient with the mgcv package in R program. Because the benthic index of biological integrity (B-IBI), generally the index scores as health assessment for each stream in this region [19] , combined several univariate community metrics into an aggregate index, the B-IBI index selected as response variables for GAM regression. Furtherly, we used the GAMs because the graphical evaluation of scatterplots revealed nonlinear patterns, GAMs are well suited for fitting nonlinear response relationships that is not known a priori [23] .
The TITAN was used to analysis the magnitude, direction and uncertainty of responses of individual taxa to each stress gradient in R program using TITAN2 package written by Baker and King [11] . Taxa abundances were log(x+1) transformed to reduce the influence of highly variable taxa on indicator score calculations in each data set. Taxa with the occurrence < 5% of the samples were deleted to remove outliers and to minimize the effect of a potential operator bias [12] . Briefly, according to the direction of the response into z− taxa with a negative response, and z+ taxa with a positive response, TITAN splits sample units into 2 groups along the disturbance gradient. Bootstrapping is used to test the purity and reliable threshold indicator taxa (proportion of bootstrap replicates with P ≤ 0.05 and P ≤ 0.05) and the uncertain locations of taxa and community change points [24] . The sum of the IndVal z-or z+ scores were then used to identity a potential predictor value representing the negative and positive community thresholds, respectively [11] .

Characteristic of water quality parameters
With reference to environmental quality standards for surface water in China (GB3838-2002), the majority of oxygen concentration was beyond 6 mg/L and better than standard II. However, organic pollution was serious for the relatively higher BOD 5 and COD Mn . The concentration of TP was lower, however, the higher TN and NH 3 -N which beyond the standard level III indicated a severer nitrogen pollution in this basin (Table 1). The PCA loading for each environmental parameter in the first four axes and the variance explained by each parameter were presented in the Table 2. The first four components explained approximately 87.9% of the accumulative total variance and the explained variance from the first to forth axis was 34.2%, 23.1%, 16.77% and 13.79%. Salinity, nutrient and organic pollutions were the most significant parameters loaded in the first four PCA axes. The first component was strong positive to EC, TN and NH3-N. The second component loaded with BOD 5 and COD Mn showed a strong positive relation.

Effects of water quality parameters on macroinvertebrate community
A total of 205 taxa were identified and belonged to 10 classes, 30 orders and 131 families. Throughout the whole river basin, macroinvertebrate community was dominated by Insecta (81.9% of the richness), followed by Gastropoda (8.0%), Hirudinea (3.5%) and Oligochaeta (2.6%). The results of CCA showed that the first four axes explained 60.4 % of the total variance (Table 3). Based on the forward selection criterion, CODMn, EC, and NH 3 -N were selected as the strongest parameters with the explanation ability for the variance.

Estimating water quality thresholds for macroinvertebrate community
The results of GAM showed a significantly nonlinear response of B-IBI to the axes scores of PCA, the first two axes scores could explain 78% and 66% of the variation of B-IBI, respectively (Figure 1).

Estimating land use thresholds for macroinvertebrate community
GAM showed a large change in B-IBI over the range of CropArea >25~35%, the ImperviousArea >0~5% and VegetationArea >40 ~ 70%, and could explain 80% and 81% of the variation of B-IBI, respectively ( Figure 1). The results of TITAN indicated that 53 of 68 taxa showed significant indicator value for CropArea and 50 taxa presented significantly decline with the increased CropArea%, and the taxa decline sharply among 23%~40% of CropArea (Figure 2). The change point of negative sum (z-) score was at 25.82% of CropArea (90% CI 23.38 and 33.94) (Table 4). Therefore, the declines of the majority of taxa were consistent with a community threshold. Notably, three taxa named Bithynia sp., Oligochaeta sp., Branchiura sowerbyi, increased in occurrence and relative abundance along the CropArea at 33.35% (90% CI 15.83 and 75.39).  The results of TITAN indicated that 47 of 68 taxa showed significant indicator value for impervious areas and 43 taxa presented significantly decline with the increased % ImperviousArea, and the majority of taxa decline sharply among 0.8%~7.5% of ImperviousArea. The change point of negative sum (z-) score was at 3.06% of ImperviousArea (90% CI 1.95 and 8.82). Similarly, four taxa increased in occurrence and relative abundance along the impervious areas' gradient at 2.76% of ImperviousArea (90% CI 1.32 and 19.83).

ICESD 2021
The results TITAN presented that 57 of the total 68 taxa showed significant indicator value for vegetation land. However, only identified significant declines in 3 taxa associated with vegetation land, and these three taxa increased significantly along the CropArea land gradient coincidently. The occurrence and relative abundance of these three taxa significantly decrease along the VegetationArea at 61.30% (90% CI 46.70 and 81.62). In contrast, 54 taxa presented significantly increased with increasing VegetationArea, and the taxa incline sharply among 62%-85% of VegetationArea. The change point of positive sum (z-) score was at 65.76% of VegetationArea (90% CI 61.86 and 74.16).

Discussion
Although biocenosis of stream ecosystem primarily depends on water quality , an upper limits of macroinvertebrate diversity responded to water quality were lack of enough evidence. In current study, GAM indicated that these water quality parameters limited assemblage integrity reaching at EC=162.3~422.2 μS/cm; NH 3 -N=0.19~0.33 mg/L; COD Mn =2.20~4.41 mg/L respectively. TITAN revealed significant declines in occurrence and abundance of intolerant macroinvertebrate taxa in response to rising these water quality parameters (EC=269.3 μS/cm; NH 3 -N=0.32 mg/L; COD Mn =2.40 mg/L respectively). High EC (283 μS/cm) have been proved to limit the diversity of intolerant macroinvertebrate taxa [25] , which was consistent with U.S. EPA also recommended the benchmark value of conductivity should be set at 300 μS/cm to protect the stream biota [26] . Similarly, high nutrient concentrations (NH 3 -N>0.12 mg/L) had also been considered limiting macroinvertebrate diversity in streams [12] , and ecological status would be limited when NH 3 -N beyond 0.03 mg/L . Although our results were higher compared to the results of previous studies, the GAM was similar to the change points for the negative response identified by TITAN, and controlled the NH 3 -N<0.32 mg/L that could protect majority biota. Furthermore, the threshold of NH 3 -N was superior to environmental quality standards (GB3838-2002) for II Class (<0.5 mg/L). Our results also found that the assemblage integrity and intolerant macroinvertebrate taxa would be limited when the COD Mn reached beyond 2.40 mg/L. As a complex water quality parameter, COD Mn refers to organic pollutants, which would be difficult to control its content like other single pollutant. Moreover, there was also lack of evidence for its threshold identification, yet the results based on our analysis were superior II Class (<4 mg/L, GB3838-2002).
CropArea, impervious land use and natural vegetation land loss were limiting the ecological status of macroinvertebrate community integrity, and a large change trend was found over a range of the CropArea >25~35%, the ImperviousAreas >0~5% and the loss of VegetationArea >30~60%. Threshold values reported in the previous studies showed that biodiversity declines abruptly when native vegetation loss exceeds 60% of original level [2] , and reduced across the 88.2% of the biome where less than 30% forest cover remains [27]. Although these studies were based on tropical rainforest ecosystem, disruptions involving linkages between terrestrial and stream habitats that aquatic and terrestrial organisms were affected simultaneously when the habitat loss exceeded the threshold, thus there was a roughly common threshold for both freshwater and terrestrial environments [6] .
Previous studies suggested that a similarity of responses among different groups could vary between 25 and 50% in environments under pressure from CropArea [28,29] . And numerous taxa (z-) significantly declined between 0.5 and 2% ImperviousArea, whereas change points for taxa (z+) were distributed from 1 to 33% [24] . Kail et al. also suggested that the threshold for change points of taxon-specific was 4.0% (z-) and 7.1% (z+) ImperviousArea, respectively [12] . Compared with TITAN, GAMs which derived for aggregated response variables (metrics) was similar or slightly higher. King and Baker indicated that aggregated response metrics used in biomonitoring community change were relatively insensitive in nonlinear response to environmental gradients, and aggregated community metric failed to distinguish the synchronous sharp decline of sensitive taxa because species sensitive to anthropogenic disturbance were often relatively uncommon [24] . Eventually, the threshold based on aggregated response metrics would be lagging behind taxon-specific results. However, the aggregated community metric would be considered as a threshold for ecosystem functioning and the thresholds derived from taxon-specific as change points in species conservation [12] , combining different data analysis approaches that separate negative and positive responses within assemblages will help resource managers to identify and estimate dominant stress responders within the assemblage.