The response of benthic diatom community to anthropogenic eutrophication of a river basin under agricultural influence in NE China

Abstract Anthropogenic eutrophication has universally threatened river ecological health and has been a key issue in river conservation. Agricultural sewage, which leads to increased nutrient levels and results in the loss of ecological function and biodiversity, results from the growth of agricultural activity and human populations. Consequently, attention has increasingly been focused on developing a relevant tool to describe the relationship between the biotic communities and anthropogenic eutrophication. The benthic diatom community has been widely used as an effective tool for indicating the ecological status of river systems especially in the context of recent international Water Framework Directive policies (WFD). In this study, a statistical method consisting of trophic status index (TSI), Specific Polluosensitivity Index (IPS), Descy Index (DES), redundancy analysis (RDA), and the geographic information system (GIS) technique was conducted to explore the response of benthic diatom community to anthropogenic eutrophication. Initially, the TSI, IPS and DES were determined to differentiate the trophic gradient along a gradient in agricultural intensity. Subsequently, RDA was used to identify spatial distribution patterns of environmental parameters and benthic diatom communities. Finally, the RDA scores as calculated and spatially interpolated were applied to GIS technology to reveal the response of benthic diatom communities to anthropogenic impacts along the trophic gradient. Our results revealed that the studied basin mostly exceeded the eutrophic level. Two different patterns of diatom communities response to ecological gradient were identified based on RDA, one representing the agriculture eutrophication discharges in electrical conductivity (EC) and the other representing organic pollution as chemical oxygen demand (COD). The developed method in this study highlight that the EC was more effective than the nutrient indices in determining diatom distribution in a eutrophic agriculturally-influenced river system.


Introduction
Rivers are considered to be the major recipients of organic pollution from agricultural non-point sources (Landis 2017;Ratmaya et al. 2019;Calizza et al. 2020). Agricultural pollution have become serious threat in inland and coastal regions, with the rapid growth of population and economy (Landis 2017;Chen, Teng, et al. 2019;Jiang et al. 2019). Agricultural non-point source pollution is a global concern leading to eutrophication in rivers, which subsequently make the polluted rivers to loose their ecological function, biodiversity and to represent negative impacts on human health (Yan et al. 2016;Li et al. 2019;Calizza et al. 2020). Previous studies showed that the growing pressures on agricultural practices have led to many issues such as scarcity of water supply, chemical contamination, and salinization in temperate rivers (Chen, Teng, et al. 2019). Once agricultural influenced river systems become fully eutrophic and contaminated, it is difficult to restore to their original state. Heilongjiang Province is the most agriculturally developed region in China (Chen, Teng, et al. 2019). The Lalin River Basin (LLRB), situated in the south part of Heilongjiang, covers an area of 19923 km 2 and annual runoff is 350 million m 3 . The LLRB is a major drainage basin of the Songhua River with nutrient levels significantly influenced by agriculture practices (Yin et al. 2018). More than half of the basin is occupied by rice and corn cultivation. The rapid development of rice and corn cultivation in LLRB has resulted in the rapid deterioration of ecological function and biodiversity in the LLRB (Yin et al. 2018;Chen, Teng, et al. 2019). Previous studies have shown a significant increase in water nutrition and pollution along the gradient of agricultural activity intensity in LLRB (Liu et al. 2013;Yin et al. 2018).
Phytoplankton is widely distributed in marine and freshwater ecosystems (Chen et al. 2016;Jiang et al. 2019). Biological and ecological traits of phytoplankton communities are considered relevant tools to indicate environmental conditions in the water column (Chen et al. 2016;Gang et al. 2019;McQuatters-Gollop et al. 2019). The benthic diatom community has ecological preferences and biological attributes filtered by selected adaptations and habitat constraints (Gudmundsdottir et al. 2013;Chen et al. 2016;Chen, Zhang, et al. 2019). The distribution pattern of benthic diatom communities is influenced by organic pollution, eutrophication, flow velocity, and other factors (Mangadze et al. 2017;Chen, Zhang, et al. 2019). They are useful indicators because of their short life-history and quick response to environmental changes (Schneider et al. 2012;Kafouris et al. 2019). Recent studies showed that diatoms communities and indicator indices can indicate many environmental influences such as acidification, heavy metal pollution, land use and eutrophication (Bere et al. 2016;Yan et al. 2016;Petrou et al. 2019;Pillsbury et al. 2019). Moreover, WDF (Water Framework Directive, WFD) using diatom communities to assess the ecological status, has shown good results (Karaouzas et al. 2019). However, the pattern of benthic diatom communities, indicate indices, and their relationship to environmental parameters in LLRB have been poorly investigated, although previous studies in other rivers from similar latitude have highlighted the value diatom indicators for some anthropogenic pressures (Becker et al. 2018;Pardo et al. 2018;Karaouzas et al. 2019) Development of effective approaches to understand the relationship of biotic communities and agricultural eutrophication is crucial to better understand and control the agricultural pollution in aquatic ecosystems (Wang et al. 2018). For this purpose, several statistical approaches such as unconstrained ordination, Principle component analysis (PCA) and Nonmetric multidimensional scaling (nNMS) have been widely conducted (Gu and Gao 2019;Petrou et al. 2019). Self-organizing Maps (SOM), while Cluster analysis (CA) were performed to explore the plankton pattern response to anthropogenic eutrophication (Lu et al. 2017;Wang et al. 2017). Trophic State Index (TSI) is a classification system designed to 'rate' individual lakes, ponds and reservoirs based on the amount of biological productivity occurring in the water (Ludovisi and Poletti 2003;Adamovich et al. 2016). Redundancy analysis (RDA) is a typical constrained ordination approach containing species and environmental weighted data and is frequently used in aquatic ecological surveys (Chen et al. 2016;Jiang et al. 2019). Spatially interpolated data is a useful method to reveal the gradient of anthropogenic pollution in the water column and sediments (Gu and Gao 2019). The spatial distribution of RDA weight is useful to reveal the relationship of biotic communities to environmental patterns (Pillsbury et al. 2019). In this context, a method, based on TSI, RDA and GIS, was developed to explore the response of benthic diatom communities to anthropogenic eutrophication along an agricultural practice gradient in the LLRB. Namely, TSI was applied to determine the degree of trophic status by nutrient input from agricultural practice in LLRB. Subsequently, RDA coupled with GIS, was performed to export continuous spatial patterns at regional scales and to reveal the response of diatom communities to anthropogenic eutrophication. We hypothesized that the spatial distribution pattern of benthic diatom communities would closely couple with the trophic level in agricultural influenced river ecosystem. When the spatial distribution pattern of dominate species, indicator indices and Shannon-weaver index couple with the TSI and environmental factors, we can verify our hypothesize. The aims of this study were: (1) to investigate benthic diatom community in a typical agricultural influenced basin in China firstly; and (2) to explore the effects of environmental factors on benthic diatom community distribution pattern based on multivariate statistical analysis.

Study area and sampling
Since the 1970s, growing urbanization, population, industries and agricultural practice have severely impacted the rivers, lakes and wetlands of northeast China (Chen, Teng, et al. 2019). LLRB, covering an area 19923 km 2 , is located east of the Songhua River (125 34 0 $128 34 0 E and 44 00 0 $45 30 0 N). The mean annual rainfall of $500-800 mm. LLRB is famous for its huge grain production. Midstream and downstream areas of the LLRB are mainly cropland within an area of about 8400 km 2 . Annual discharge of COD, NH 4 þ , TN and TP generated by domestic sewage, garbage, livestock manure, chemical fertilizer and pesticide in LLRB is 100000 tons, 16900 tons, 29600 tons and 2881 tons respectively (Yin et al. 2018). Former studies showed that the impact of anthropogenic activities on water and soil in LLRB have severely impacted the river (Liu et al. 2013;Chen, Teng, et al. 2019). Concentrations of Co and Se were slightly higher than their background values, meanwhile, Cd and Ni in some areas of LLRB were greater than the Grade II of the Standard according to the Soil Environmental Quality of China (Chen, Teng, et al. 2019). The midstream water quality in downstream sections of LLRB ranged from class III to V according to the Water Environmental Quality of China. Altogether, as one of the most densely rice and corn cultivated regions in China, aquatic environment has been severely impacted by anthropogenic activities such as agriculture sewage discharge (Chen, Teng, et al. 2019).
For this study, benthic diatom and water samples were collected in rainy season (august) 2018 from twenty-nine sampling stations distributed throughout LLRB ( Figure 1). Benthic diatom samples were preserved with 1% Lugol's solution (Chen et al. 2016). The diatom samples were stored at 4 C for 24 hours until treatment with 30% H 2 O 2 to remove organic material and mounting on slides with Naphrax (Potapova and Charles, 2007). Diatom counts from each station were expressed as relative abundances to eliminate the effects of rare species. At least 600 diatom valves were counted and identified employing a Zessis Imager A2 microscope at 1000 Â magnification (Chen, Zhang, et al. 2019). All diatom taxa with less than two occurrences and a percentage < 5% were removed from database (Karaouzas et al. 2019).
Species with a relative abundance of more than 5% in at least one sample were determined as the dominant species in this study (Chen, Zhang, et al. 2019). A total of 9 environmental parameters were measured in this study. Physical parameters such as water temperature (WT), Dissolved oxygen (DO), and pH were conducted by YSI (Perkin Elmer Inc., USA) in situ. Nutrient analyses included COD, TN and TP and chlorophyll-a concentration were measured according to Chinese national standards for water quality (Wang 2002; similar to those of the American Public Health Association (Gilcreas 1966)).

Multivariate statistical analysis
In order to explore the response of diatom communities to environmental change and determine the key factors, constrained ordination and GIS-based approach were conducted between 21 dominant diatom taxa and 9 environmental parameters. The initial analysis was to determine the comprehensive trophic state index (TSI) to distinguish the ecological gradient such as trophic levels. In the second analysis, RDA was performed on 21 dominant diatom species and 9 environmental parameters to explore the response of diatom communities to the ecological gradient. The last analysis, GIS approach (Inverse distance weight, IDW) was conducted to visualize the response of benthic diatom communities to anthropogenic impacts based on RDA scores. According to Leps (Leps et al. 2003) the DCA analysis showed that the gradient lengths were less than 2 standard deviation units, indicating that a linear model (e.g., RDA) can be used (Leps et al. 2003). Monte Carlo permutation test was performed to select the significant variables (499 random). IDW approach was carried out to illuminate spatial pattern of RDA scores by geo-statistical analysis module of GIS (Gu and Gao 2019). Prior to multivariate statistical analysis all variables including species and environmental parameter data were subjected to a log (X þ 1) formula to evaluate the normality.
The equations for TSI are calculated as follows:  (Wang 2002).
To assess the trophic levels, Specific Polluosensitivity Index (IPS) and Descy Index (DES) were determined. The IPS and DES were performed by Package 'diathor' by R (version 4.2.1). The lower value of IPS and DES indicating the higher trophic levels.
One-way analysis of variance (ANOVA) was performed to test for determined differences in 9 environmental parameters among different groups (95% confidence intervals). Pearson correlation coefficients was performed to test the correlation between Shannonweaver index and 9 environmental parameters. Analysis of Similarities (ANOSIM) and Similarity Percentages-species contributions (SIMPER) were conducted to explore the difference of diatom communities and filter the main dissimilarity contributed species. ANOVA and Pearson correlation coefficients were determined using SPSS 22.0 software; ANOSIM and SIMPER were performed by PRIMER 5.0 (Gibert and Escarguel 2019;da Silva et al. 2020). RDA and DCA were performed by Canoco for windows 4.5 software. IDW were carried out by ArcMap 10.3. The land redundant use information such as cropland was downloaded from remote sensing data (Resource and Environmental Data Cloud Platfsorm http://www.resdc.cn/).

Environmental parameters spatial characteristics
All nine environmental parameters displayed wide variations except WT (Figure 2). There was a clear rising tendency in TSI from upstream to downstream (Figure 3). TSI indicated that sampling sites were mainly eutrophic except station 16 in the study period. Combined remote sensing data and TSI found that three different groups of all 29 stations were identified in LLRB according to their trophic level and agriculture practice (Figure 1b). Group-1 included station 1 to station 4, station 16 and station 23, which were characteristic by the meso-trophic to light eutrophic. Group-2 included stations 5 to 12, station 17 and station 25, which were characterized as meso-trophic. Stations 13 to 15 and stations 26 to 29 belonged to Group-3, which was characterized as highly eutrophic. ANOVA showed that there was substantial difference of EC, COD, TN and TP among three groups (p < 0.05) (Table 1). DO, pH, EC, ORP, COD, TN and TP were found to have similar distribution characteristic and their highest values were mainly recorded at Group-2 and Group-3. Contrary the highest value of turbidity was recorded at Group-1.

Benthic diatom communities and indices
A total of 130 taxa belonging to 30 genera were identified in the LLRB. Five taxa belong to Centricae, and 125 taxa belong to Pennatae. The highest species richness was mainly observed in Navicula Bory. (16), Gomphonema Ehrenberg (16), Nitzschia Hassall (15), Surirella Turpin (9) and Fragilaria Lyngbye (8). 21 dominant diatom taxa including Melosira varians and Encyonema minutum were identified (Figure 4). Species richness ranged from 11 in station 8 to 35 in station 28. ANOVA indicated significant differences among three groups for Chlorophyll-a (p < 0.05). The higher Chlorophyll-a concentrations were mainly recorded in group-3 while the lower Chlorophyll-a concentration were commonly recorded in group-1 (Figure 4). Richness of groups 1, 2, 3 was 59, 51 and 50 respectively. Shannon-weaver index at sampling stations along group-1 to group-3 was slightly higher, with means of 2.09, 2.24 and 2.55 respectively (Figure 4). Twenty one dominant taxa were identified. Analysis of similarities (ANOSIM) revealed significant difference between group-1 and group-3 in the diatoms communities (Global test R ¼ 0.13, p < 0.05). Similarity percentages-species contributions (SIMPER) indicated that average dissimilarity between group-1 and group-3 was 68.53 (Table 2). Rhe results of SIMPER indicating that Melosira varians (23.64%) Gyrosigma scalproides (13.02%), Navicula radiosiola (9.48%) and Nitzschia palea (6.42%) were main dissimilarity contributing taxa of three groups ( Table 2). The spatial distribution pattern of IPS and DES were similarity (Table 3). ANOVA showed there was significant different among three groups. Mean lowest value of IPS and DES were recorded in Group-3, indicating the higher trophic level in this region. The results of IPS and DES highlight the three groups, which divided by TSI could indicating the trophic gradient in LLRB.

Multivariate statistical analysis
The results of RDA showed that the first and second ordination axes explained 37.4% and 18.5% of the total variance, respectively. A Monte Carlo permutation test indicated that EC and COD were the key factors driving benthic diatom communities ecological distribution. Axis 1 was positively correlated with EC (0.8843). Axis 2 was negatively correlated with COD (0.6573). A higher level of EC was associated with higher Navicula cryptotenella, Halamphora montana, and Gyrosigma scalproides. A higher level of COD was associated with higher Navicula radiosiola, Nitzschia palea and Gomphonema parvulum var. lagenula. RDA sample score was visualized with IDW technology by GIS approach (Figure 5). In current study, there was no significant correlation between Shannon-weaver index and environmental parameters (p > 0.05) ( Figure 6). Shanon-weaver index was negatively correlated with Tur.(p > 0.05) and DO, and positive correlated with another 7 environmental parameters (p > 0.05).

Discussion
In recent years, nitrogen and phosphorus leaching from agricultural practice lead to river eutrophication and pollution had been recorded extensively in developing countries Ratmaya et al. 2019). China has consistently become one of the largest producers of the world's organic agricultural products over the past decade. The Lalin River Basin is a major tributary of the Songhua River of northern China, which is significantly influenced by anthropogenic activities such as agriculture (Yin et al. 2018).

Benthic diatom communities as an indicator in eutrophic agriculture river
In current study, we provide the key information about the he distribution of benthic diatom communities including, taxa richness, relative abundance, Shannon-weaver index, and dominant taxa in LLRB firstly. A total of 130 diatom taxa were identified, the diatom   communities were similar to same latitude Spain rivers (132 taxa, 676 sites) (Pardo et al. 2018). In general, Shannon-weaver index was widely used to reflect the community heterogeneity in local scale. In this study, the Shannon-weaver index of diatom communities ranged from 1.13 to 3.12 means 2.29, suggested a lower ability to tolerate anthropogenic eutrophication stress in LLRB. The correlation analysis showed that Shannon-weaver index was negatively relate to TP, indicating that eutrophic reduced the diatom community stability. Based on TSI all 29 sampling stations were divided into three categories. The meso-trophic to highly eutrophic nutrient states that can provide the greatest competitive microhabitats for eutrophic and pollution tolerant diatom taxon, particularly Melosira, Navicula, Encyonema and Nitzschia species (Chen et al. 2016;Wang et al. 2017;). By combing the information from TSI, IPS, DES and ANOSIM, we found that diatom communities in Group-1 were significantly different than Group-3, although all groups were mainly dominated by pollution or eutrophication tolerant diatom taxa. These results may imply that benthic diatoms are significant signals in response to environmental filters in eutrophic river system (Wood et al. 2019). A clear separation of diatom taxa was observed in the RDA ordination indicating that distribution of benthic diatoms is closely related to ecological gradients along nutrient levels ( Figure 5-a). In group-1, which is characterized as meso-trophic to eutrophic, the diatom communities were dominated by Melosira varians and Encyonema minutum (Figures 4 and 5-a). Melosira varians which was dominant in group-1 was considered a good indicator of eutrophic habitats that occur in lowland streams characteristic of meso-trophic conditions in South Africa and affected by agriculture (Petersen et al. 2018). Encyonema is widely accepted as a freshwater genus frequently associated with eutrophic or nourishing conditions. For example, Encyonema minutum has been described as a polluted water indicator (Lavoie et al. 2008). Prior study recorded that Encyonema minutum was closely related to the higher nutrient concentrations origin from agriculture and sewage water in Andean river basin. The high dominance of Melosira varians and Encyonema minutum in group-1 is indicative of eutrophic status here. Our TSI also proved that a meso-trophic to light eutrophic state in Group-1. The aquatic trophic gradient was considered to constrain on the diatom assembly directly by change the in nutrient concentrations . Moreover, using benthic diatoms to indicate the trophic gradient has been reported in numerous aquatic systems (Gudmundsdottir et al. 2013;Vilar et al. 2014;Stenger-Kov acs et al. 2018). In this study, Group-3 exhibited the highest chlorophyll-a concentration and TSI, and was dominated by Nitzschia palea, Navicula cryptotenella and Gyrosigma scalproides. Several studies have recorded that Nitzschia species densities may increase as a consequence of nutrient enrichment in river systems Pardo et al. 2018).
Prior study noted enrichment promoted the growth of Nitzschia palea in nutrient-rich lowland rivers (Licursi and Gomez 2009). Additionally, Nitzschia palea and Navicula cryptotenella were considered as eutrophication and polysaprobic taxa indicators, respectively . Based on ecological attributes of dominant taxa, we further concluded that the trophic level was higher in group-3 than in the other groups. Ecological studies focusing on Gyrosigma scalproides as an agricultural pollution indicator in river system are scare. Notablely, in this study, ANOSIM and SIMPER showed that Gyrosigma scalproides was the main contributor for the dissimilarity between group-1 and group-3. In agreement with our findings, a few Gyrosigma species are considered as halophilic taxa, capable of tolerating meso-trophic conditions (Mrozinska and Czerwik 1996).
Gyrosigma was recorded as tolerant to alkaline and higher ionic condition (Ohtsuka 2005). In Group-3, pH means was 8.87 (Table 1). Meanwhile, the EC of Groups 3 was significant higher than Group1 and Group2 (Table 1). Gyrosigma scalproides in Group-3 is probably due to its ability to tolerate elevated alkaline and EC conditions. By combining with TSI and Chl-a concentration in Group-3, which exhibited the highest nutrient concentration, we infer that Gyrosigma scalproides is a potential indicator to indicate changes in nutrient concentrations along trophic levels in temperate rivers influenced by anthropogenic eutrophication. However, the ecological information of Gyrosigma scalproides on eutrophic rivers is insufficient, thus, a further survey is necessary.

Impacts of environmental parameter on diatom communities
Based on TSI, IPS and DES a major ecological gradient were determined that trophic level was observed significant change from group-1 to group-3. RDA indicated that EC and COD were important explanatory variables associated with the distribution of diatom communities in LLRB. The RDA-1 was associated with EC (r 2 ¼0.567, p < 0.05), representing the pattern of diatom communities along EC gradient; RDA-2 was associated with COD, likely representing organic pollution impacted on diatom communities. Both of two patterns (RDA-1 and RDA-2) of diatom matrix response to anthropogenic impacts were visualized by IDW in Figure 5. Previous study note that EC is an relevant indicator to indicate the total ionic concentration in aquatic system (Ryves et al. 2002), and is a substantial factor to dynamic diatom communities distribution in wetlands (Tibby et al. 2007), peathlands (Chen et al. 2016) and Rivers . Former studies have reported EC strongly impact the diatom community assembly in many regions and similarly had a significant ecological influence in our study (Potapova and Charles 2003).
Other evidence for the influence of EC on diatom communities was reported in the Bloukrans River, researcher founded that diatom communities is reliable bio-indicator for assessment conductivity in small austral temperate river system (Mangadze et al. 2017). For instance, Guadiana River stations dominated by Nitzschia sp. taxa including N. capitellata, N. umbonata and N. palea were influenced by by EC caused by agriculture practice (Urrea and Sabater 2009). Additionally, Pestryakova et al. (2018;Pestryakova et al. 2018) reported that high relative abundances of Melosira, Navicula and Nitzschia taxa, which were also found in high EC in Yakutian Lake. In our study, EC varied from 59.67 lS/cm to 404.47 lS/cm, is consistent with Yakutian lakes (Pestryakova et al. 2018) and Bloukrans River (Mangadze et al. 2017) suggest water extremely vigorous in electrolytes, similarly to some polluted river in America (Potapova and Charles 2003). Unfortunately, no information of total ionic is available. As expected, the 'hot spot' of RDA-1 showed that the pattern of diatom communities influenced by EC were conspicuous in group-3 ( Figure 5-b). According to the Report on the State of the Ecology and Environment in China (http://english.mee.gov.cn/Resources/Reports/soe/) in the Songhua river basin in 2018 and cropland distribution data from remote sensing, the higher EC of LLRB are mainly located in group-3, and associated with agriculture eutrophication discharges . The consistency of Navicula, and Nitzschia taxa affected by EC in LLRB highlights that EC was a substantial factor dynamic diatom communities in agriculture river system. COD was considered as an organic pollution signal in freshwater environmental monitoring, meanwhile, is severely influential to diatom community assembly in streams (Castillejo et al. 2018). Organic pollutants that are discharged into rivers may through oxidation expand the valve of diatoms, that cause diatom teratologies (Lavoie et al. 2017). Recent studies of bio-assessment using diatom assemblage in rivers have also proved the importance of organic pollution in benthic diatom community assembly (Bere and Tundisi 2011;Wood et al. 2019). In our study, IDW showed that all though the RDA-2 'hot spot' was part located in group-3 but the spot were mainly separation, this result suggested that the response of diatom communities to COD gradient is relative fuzzy (Figure 5-c). TN and TP were widely used as dynamic factors to determine the plankton and microorganism community assembly. In this study, ANOVA indicated that the other nutrient index such as TN and TP were significant different among three groups. But Monte Carlo Permutation test revealed that the correlation of TN and TP with diatom communities is not significant (p > 0.05). Combined with the separation IDW 'hot spot' (RDA-2), we implied that in eutrophic agriculture river system the degree of explanation by nutrient indices maybe weaker. This result may be attributed the higher background nutrient level. On the contrary, based on information from TSI, RDA and GIS approach we found that the EC is more relevant indicator to dynamic diatom communities in eutrophic agriculture river system. In our count, we found some teratologies that were distributed regularly along trophic level, such as Navicula cryptotenella, Gyrosigma scalproides and Gomphonema parvulum. Unfortunately, in this study we didn't test the response of teratologies of diatoms communities to trophic level. In future, the develop the study focus on the succession of diatom community between cropping season and ley farming season is benefit for the water quality assessment in agricultural practice river ecosystem.

Conclusion
Our results revealed that different nutrient state condition of agricultural influenced rivers could be distinguished by diatom communities. EC and COD were significant environmental factors associated with diatom community assembly. The interpolation analysis based on GIS in this study suggested EC had a better effective to describe the diatom spatial distribution. In view of the visualization and interpolation of GIS approach, the method include TSI, IPS, DES, constrained ordination and GIS used for routine rivers ecological study is worth to further explored. Overall, our results highlight the utility of diatom communities in the environmental assessment of agricultural influenced basin and developed a comprehensive approach to reveal that diatoms response to environmental gradient.