Using eDNA to Identify the Dynamic Evolution of Multi-Trophic Communities Under the Eco-Hydrological Changes in River

As significant players in material cycling and energy flow, bacteria and eukaryotes play a vital role in the ecosystem. Nevertheless, the community dynamics of bacteria and eukaryotes in rivers and their responses to changes in ecological hydrology have not been studied thoroughly. Based on eDNA technology, this study investigated the bacterial and eukaryotic communities in the upper, middle and lower reaches of the Weihe River in different seasons. The seasonal variation and geographical distribution of bacterial and eukaryotic community structures showed significant heterogeneity. The selective theory well explained the response of microbial community assembly to seasonal changes. Deterministic processes dominate microbial community assembly in the middle and lower reaches. The composition and metabolic potential of key functional genes of nitrogen and phosphorus cycling (nosZ, pqqB, pqqD, and pqqE) exhibited strong seasonal patterns and were significantly correlated with the physical and chemical properties of water. There were significant differences in molecular ecological networks in different periods (p < 0.05), with a gradually increasing trend in the complexity of the network from winter to summer. The keystone species (Hub) of the microbial food web in each season included microorganisms (Malikia), algae (Stephanodiscus), and invertebrates (Polyarthra). Structural equation modeling (SEM) results indicated that invertebrate was an important driving factor affecting the changes in community structures. In micro-food webs, both “bottom-up” (resources) and “top-down” (predation) forces strictly controlled the relationship between taxa. Nitrogen and phosphorus concentrations affected microbial networks, and there was a significant correlation between bacterial and eukaryotic groups and eco-hydrological variables (p < 0.05). Furthermore, we identified the taxon’s change point using threshold indicator taxa analysis (TITAN), quantitatively revealing the response thresholds of taxa to eco-hydrological changes.


INTRODUCTION
Climate change and human activities significantly reduce the river's resilience to disturbance, and aquatic biodiversity controls the basic ecological processes and functions of the river . As major players in the material cycle and energy flow, bacteria and eukaryotic communities are crucial in driving ecosystem processes . Dynamic changes in microbial communities exert direct and profound impacts on the ecological integrity of freshwater ecosystems (Tang et al., 2020). Such changes are deemed essential to understanding the rapid changes in the structure and function of aquatic ecosystems (Niu et al., 2018). Researchers have revealed how human activity affects microbes' functional changes and metabolic potential by describing the spatial distribution patterns of microbial and eukaryotic, including fluctuations in diversity, dispersal patterns, and levels of species interactions in different environments (Webster et al., 2018). Related studies showed a close correlation between bacteria and eukaryotic community composition (Logares et al., 2018). Their diversity was influenced by land-use practices, environmental factors, hydrological properties, and geographic distance between sampling sites (Xie et al., 2017). However, challenges remain in studying the interactions of a wide range of bacterial and eukaryotic taxa across different trophic levels in rivers with strong anthropogenic impacts (Hu et al., 2022).
Environmental DNA (eDNA) can extract the community composition information from the samples by expanding the target DNA fragments, suitable for environments with rich biodiversity, such as rivers (Suter et al., 2021). eDNA technology can simultaneously investigate organisms at multiple trophic levels, providing critical information about the complex biological interactions associated with ecosystem changes . The eDNA technology has allowed for accurate monitoring of organisms, effectively improving the timeliness of biological evaluation and validating the potential of the technology to assess the status of aquatic ecosystems (Bagley et al., 2019). The 16S rRNA and 18S rRNA genes have high species coverage and relatively complete sequence databases (Walters et al., 2016;Xu et al., 2020;Sun et al., 2021). These gene markers have been widely applied to monitor the species composition and community structure of marine and freshwater bacteria and eukaryotes (Harper et al., 2019). In addition, some researchers have studied the seasonal changes and spatial composition of eukaryotic phytoplankton diversity based on the 18SrRNA gene and further revealed the impact of environmental variables on phytoplankton functional groups and community composition (Yang et al., 2021a).
However, these studies may fail to accurately describe the response of aquatic communities to environmental changes. To explore the mostly non-linear relationships between stress factors and biological parameters, ecologists have proposed non-linear thresholds analysis methods, such as nonparametric changepoint analysis (nCPA), threshold indicator taxa analysis (TITAN), and random forest (Sultana et al., 2019;Triadó-Margarit et al., 2019). These non-linear methods show higher accuracy in building models to resolve biological-environmental data relationships than the linear models. Especially, TITAN can detect the critical point where the distribution law of each community's characteristic indicators changes along an environmental gradient (Martins et al., 2021). Further, it can assess the synchrony of the crucial issue of different community characteristic indicators and determine the critical point where the whole community responds to environmental factors. Recently, the application of TITAN to identify thresholds with habitat changes and habitat conditions improved the ability to balance resource conservation and development while identifying sensitive indicators (Xuan et al., 2019).
Revealing the dynamic evolution mechanism of communities at different trophic levels is of great scientific significance and practical urgency to guarantee water safety and ecological security in river basins. At present, the ecological environment of the Yellow River basin is fragile, the water security situation is complex, and the quality of development needs to be improved (Yang et al., 2021b). As the largest tributary of the Yellow River, the Weihe River also faces problems such as a sharp decline in biodiversity, an increase in hydrological rate, and severe water pollution (Zhao et al., 2020). However, existing studies in the Yellow River Basin have only focused either on revealing individual organismal groups or single functions, seriously ignoring the in-depth understanding of the structures and functions of the multi-trophic communities in the river under changing environments .
In addition, they have yet to quantify the response relationship between aquatic communities and eco-hydrological changes . Therefore, based on eDNA metabarcoding technology, we explored the dynamic evolution mechanism of biodiversity from bacteria to invertebrates. This study aims to 1) characterize the spatiotemporal changes in the structure and function of bacterial and eukaryotic communities in Weihe River, 2) elucidate the coexistence phenomenon and potential interactions between coexisting taxa at different trophic levels, and 3) quantitatively analyze the threshold of biocommunity response to eco-hydrological variables.

Study Area
The Weihe River basin covers an area of 134,766 km 2 with complex geomorphological features, spanning three different landform units: the Loess Plateau, the Guanzhong Basin, and the Qinling Mountains (Song et al., 2018). The annual rainfall of the Weihe River is between 350 and 700 mm, concentrated in summer, with an apparent seasonal difference (Song et al., 2021). In recent years, due to the rapid development of the Guanzhong urban agglomeration, the water pollution in the Weihe River Basin has intensified, the biodiversity declining sharply, and the ecological function of the river deteriorating (Liang et al., 2022).
When sampling bacteria and eukaryotes, the operation person wore masks and disposable medical gloves to collect 2,000 ml of water samples from each water layer of the monitoring sections and cryopreserved them until transported to the laboratory (Bahram et al., 2021). Under sterile conditions, a 0.45 μm micropore filter was applied to filter eukaryotes, and a membrane filter with 0.22 μm pore size to filter microorganisms from water samples (Chi et al., 2021). The filter membranes were placed in sterile centrifuge tubes and stored in a -80°C cryogenic freezer for the biome sequencing analysis.

Water Quality and Hydrology
A portable multifunctional water quality analyzer (HACH HQ40days) was used to measure water temperature (WT), pH, and dissolved oxygen (DO) on site (Wu et al., 2019). Water quality variables, i.e., total phosphorus (TP), total nitrogen (TN), and ammonia nitrogen (NH 4 + -N), were measured by spectrophotometry (State Environmental Protection Administration, 2002). Flow velocity (V), river discharge (Q), and water level (D) data were obtained by using a portable flow meter (MGG/KL-DCB). The sections located in the national hydrological and water quality monitoring stations Beidao (BD), Xianyang (XY), and Huaxian (HX) hydrological stations) adopt the real-time monitoring data of the national monitoring station.

DNA Extraction and Sequencing
According to the instruction, DNA was extracted via the Fast DNA SPIN extraction kit (MP Biomedicals, Santa Ana, CA, United States), and the samples were stored at -80°C for further analysis. The quantity and quality of the extracted DNA were measured by using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and agarose gel electrophoresis, respectively. A detailed description of the amplification protocol has been published (Valentini et al., 2016). This present study performed samples for two-end sequencing on the Illumina MiSeq PE250 platform (Illumina, San Diego, CA).

Bioinformatics and Statistical Analysis
After sequencing, the species composition, diversity index, and principal coordinates analysis (PCoA) of the microorganisms were performed using the R (version 3.5.2) (https://www.rproject.org/) "vegan" and the "ggplot2" package, QIIME (v1.8. 0) software, and the mothur software. Network analysis of the bacterial and eukaryotic communities was performed using the "Hmisc" and "igraph" packages. In order to reduce false-positive prediction and increase the sensitivity of the network, we used samples with a proportion above 0.01% for network analysis . To better understand the potential functional and metabolic pathways of specific biological taxa in water bodies, we used PICRUSt version 2.0.0 to analyze the functional content of the 16S rRNA and 18S rRNA marker genes. Compared to its original version, the reference genome database of PICRUSt 2.0.0 was more than 10-fold expanded, and its prediction results are more rigorous and accurate (Douglas et al., 2020). To analyze the direct and indirect effects of nitrogen and phosphorus nutrient changes on multi-trophic communities, we used a Structural equation modeling (SEM) to determine the direct and indirect paths among environmental change and trophic transfer efficiencies Ma et al., 2020). The SEM model was constructed using Amos (version 18.0, IBM, United States).

Identification of Threshold Responses and Indicator Taxa
Finally, TITAN was applied to determine the ecological protection thresholds of eco-hydrological indicators for species composition changes in Weihe River (Simonin et al., 2019). Before data analysis, the top 50 species samples (genus level) were selected to identify the best change point by TITAN to obtain maximum indication merit (IndVal). All species were negative response taxa (z-) and positive response taxa (z+) based on the relative abundance and frequency on both sides of the change points (Berger et al., 2018). After defining the indicator species, TITAN provides a community-level threshold reflecting the size of community change as a point of coincident change indicator of the community structure [sum (z)], indicating the presence of community thresholds if there were significant synchronous changes in the distribution of taxa (Valente-Neto et al., 2021). In addition, bootstrap technology was employed to measure the uncertainty, purity, and reliability. The thresholds and the reliability of the indicated species in various group change points were validated by the measured uncertainty (p < 0.05), purity (≥0.95), and reliability (≥0.95).

Eco-Hydrological Variables
There were significant seasonal changes in the water environment and hydrological monitoring indicators (p < 0.05) ( Table 1). As one of the primary pollutants in the Weihe River, the average concentration of TN was higher than the Class V standard limit stated in the Environmental Quality Standards for Surface Water (GB3838-2002) in all seasons. The mean flow, velocity, and water depth values were significantly higher than those in the other seasons (p < 0.05). The mean autumn and winter temperatures were much lower than those in the other seasons, while the mean autumn and winter DO and NH 4 + -N were higher than those in other seasons.

Spatial Patterns and Seasonal Changes in Bacterial and Eukaryotic Communities
A total of 2,233,453 high-quality sequences and 1,616,821 highquality sequences were retained from the 12 monitoring sections at four seasons. The average number of sequences of the 16S rRNA genes was 46,530, ranging from 32,526 to 73,569, generating a total of 74,365 ASVs. The average number of sequences of the 18S rRNA genes was 34,400, ranging from 23,192 to 49,166, generating a total of 6,680 OTUs. Major taxa in the bacterial communities were Proteobacteria, Bacteroides, Actinobacteria, Cyanobacteria, Firmicutes, Epsilonbacteraeota, Verrucomicrobia, Patescibacteria, Chloroflexi, and Acidobacteria. Meanwhile, the sequence amount of Proteobacteria was 49.94%, significantly higher than the other taxa, indicating that it was the dominant taxa. Eukaryotic communities mainly included phytoplankton algae, fungi, protozoa, invertebrate metazoans, anneldans, molluscs, etc. Algae had the most significant number of sequences and ASVs, followed by fungi and fewer protozoa and invertebrate metazoans. The relative abundance of Bacillariophyta, Ciliophora, Chlorophyta, Cryptophyta, Ascomycota, Pyrrophyta, and Rotifera was relatively high, of which Proteobacteria (29.58%) was the dominant group. Seasonal variation and geographic distribution of the bacterial and eukaryotic communities showed notable heterogeneity (p < 0.05; Figure 2). The relative abundance of Proteobacteria was higher in the middle stream and downstream than upstream during autumn and winter. The relative abundance of Proteobacteria did not vary significantly between upstream and downstream during spring and summer. The relative abundance of Bacillariophyta was low during the spring and autumn and showed a decreasing trend from upstream to downstream.
The diversity of bacterial and eukaryotic communities showed significant seasonal changes (p < 0.05; Figure 3). The PCoA analysis results showed notable differences between the bacterial and eukaryotic communities (p < 0.05; Figure 4). The dispersion degree between bacterial samples was highest in spring, while for eukaryotic community samples, it was higher in spring and summer than in autumn and winter.

Spatio-Temporal Changes of Bacterial and Eukaryotic Community Function
Microbial genes contain a large number of pollutant degradation functions and nitrogen and phosphorus metabolism functions, and the related nitrogen and phosphorus cycle function genes show strong seasonal patterns (Smith et al., 2019). Among the 36 functional genes related to the nitrogen cycle, 11 of them showed a relative abundance that was significantly affected by seasonal changes (p < 0.01; Figure 5). The relative abundance of nosZ and NORC genes used to characterize denitrification, hdh gene used to describe anammox, and napC gene used to represent dissimilatory nitrogen reduction were relatively high, showing the seasonal variation characteristics of the lowest in winter and higher in other seasons. Moreover, seasonal variation significantly affected the relative abundance of 15 among the 52 phosphorus cycle-related functional genes (p < 0.01). These 15 genes were mainly distributed in three functional groups, phosphonate mineralization, inorganic phosphate solubilization, and phosphonate transport, mainly including phnF、phnG、phnI、phnJ、phnL、phnM、phnN、phnK、 gcd、glpT、glpP、pqqC、pqqB、pqqD and pqqE genes.
A significant correlation was between the physical and chemical properties of water and the genes of nitrogen and phosphorus cycle function ( Figure 6). WT had a significant positive correlation with the denitrification gene nosZ and the dissimilatory nitrate reduction gene in nitrogen cycle function genes (p < 0.01), while it had a significant positive correlation with the inorganic phosphorus dissolution genes pqqB, pqqD, and pqqE in the phosphorus cycle genes (p < 0.01). TN and TP were significantly correlated with nitrogen and phosphorus cycle function genes (p < 0.01). The increase of TN concentration affects the expression of phosphorus cycle genes, while the rise of TP concentration affects the expression of nitrogen cycle function genes. The change of DO concentration notably affects the expression of pqqB, pqqD, and pqqE genes, indicating that dissolved oxygen would impact the biomass of some strains containing inorganic phosphorus dissolving genes. In addition, there is a significant positive correlation between nitrogen and phosphorus transformation functional genes (p < 0.001; Figure 7).

OTU Co-occurrence Network of Multi-Trophic Groups
In order to clarify which species have significant correlations, we used different colors to represent different biological groups. The species screening and significance test results indicated that the biological network comprises 592 bacteria, 40 fungi, 126 algae, and 291 invertebrates (Figure 8). A gradual increasing trend from winter to summer existed in the complexity of the network. Compared with the separately generated random network, the heterotrophic food web was more significantly clustered Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 929541 5 and forms a "small world". Further analysis of critical nodes in the bacterial and eukaryotic ecological networks revealed that nutrient distribution of taxa included decomizers, producers, and consumers.
The positive and negative correlation in the network can infer the symbiosis or competition/predation relationship between water biological communities. In this study, the positive and negative correlations of networks coexisted between different seasons, among which the more intense competition/predation phenomenon between different biological communities in the water network existed in summer. In contrast, the cooperative relationship and competitive relationship of networks coexisted in other seasons. The positive interactions in the network were mainly concentrated in modules with the same species classification, such as Cryptomonas and Chlamydomonas, showing synergy among species, while the other species were mainly negatively correlated. The species that are significantly related in multi-trophic group networks usually come from different trophic levels and a large number of edges target species with low trophic levels to those with high trophic levels or species with high to low trophic levels.
Moreover, we focused on the keystone species (Hub) of the microbial food web in each season (Figure 8), which were mainly composed of microorganisms (malikia), algae (stephanodiscus), and invertebrates (polyarthra). The abundance of bacteria was much higher than that of other groups, indicating its leading role in heterotrophic food webs.
Amoeba, which could feed on attached bacteria, was the main protozoan predator of bacteria in the river ecosystem. Arthropods constituted a high nutritional level in the microbial food web, in which rotifers were the main predators of most algae and animals in the network. The low centrality of some hub groups indicated that the network structure was less dependent on individual species.
The running results of the SEM model showed that in the micro-food webs, both "top-down" (predation) and "bottomup" (resources) forces strictly controlled the relationship   Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 929541 7 between multi-trophic groups (Figure 9). Invertebrate was a critical driving factor in the change in community structure. Heterotrophic bacteria, fungi, and algae had significant effects on invertebrates (p < 0.05), with the path coefficients of 1.01, 0.89, and 0.63, respectively, indicating that invertebrates had the most significant direct impact, followed by algae. Invertebrates had a significant (p < 0.001) effect on heterotrophic bacteria, fungi and algae with path coefficients of 1.01, 0.89 and 0.63, respectively. Algae were remarkably associated with fungi (p < 0.01), with path coefficient of 0.17. In addition, a negative correlation existed between microorganisms (heterotrophic bacteria and fungi) and the nitrogen and phosphorus cycles. Algae positively correlated with nutrients (nitrogen and phosphorus), indicating that algae promoted nutrient circulation via absorption and utilization of nitrogen and phosphorus. Biodiversity-mediated pathways had a crucial impact on ecosystem function.

Effects of Eco-hydrological Variables on Bacterial and Eukaryotic Communities
Correlation analysis ( Figure 10) showed a significant correlation between bacteria and environmental and hydrological variables (p < 0.05). Distinct positive correlation between Proteobacteria, Actinobacteria, Chloroflexi, and Acidobacteria and water temperature (p < 0.05). Hydrological variables strongly affected microorganisms, with Actinobacteria, Chloroflexi, and Acidobacteria having a notable positive correlation with the river discharge (p < 0.05). Moreover, a remarkable correlation existed between eukaryotes and environmental and hydrological variables (p < 0.05). Bacillariophyta was significantly positively associated with total phosphorus (p < 0.05), and the relative abundance of Ciliophora, Chlorophyta, Cryptophyta, Ascomycota, Rotifera, and Chordata were significantly positively associated with water temperature (p < 0.05). Bacillariophyta, Chlorophyta, Ascomycota, Pyrrophyta, and Rotifera were significantly positively associated with the river discharge (p < 0.05).

Non-Linear Responses of Bacterial and Eukaryotic Communities to Key Eco-Hydrological Variables
Based on the correlation results, the most significant two factors, water temperature, and river discharge were selected for TITAN analysis with the top 50 abundance species, considering the change points in significant indicator taxa with purity and reliability greater than 0.95 ( Figure 11, Figure 12). The TITAN results for bacterial communities indicated that the change point of WT was 24.70°C for the sum (z−) and 5.70°C for the sum (z+). The change point of river discharge was 152.79 m 3 /s for the sum (z−), whereas the change point was 12.30 m 3 /s for sum (z+), and the filtered score was 106.81 for fsum (z+). The threshold range was 18.25-24.70°C for the negative indicator taxa (Arcobacter and Flavobacterium), and 1.75-17.65°C for the positive indicator taxa (CL500-29 marine group and Mycobacterium). The threshold range for the negative response taxa (Flavobacterium and Sediminibacterium) in the runoff indication taxa was 106.80-371.79 m 3 /s, while for the positive response taxa (Arenimonas and Mycobacterium) was 8.00-106.78 m 3 /s.
The TITAN results of eukaryotic communities showed the change point of water temperature was 26.10°C for the sum (z−), and 4.25°C for the sum (z+). The change point of flow was 140.70 m 3 /s for the sum (z−), and 23.58 m 3 /s for the sum (z+). The threshold range was 8.00-29.65°Cfor the negative indicator taxa (Spumella and Stylonychia), and 0.90-8.00°C for the positive response taxa (Unruhdinium and Monoraphidium). The threshold range was 46.62-332.10 m 3 /s for the negative response taxa (Spumella and Glaucoma) in the runoff indication taxa, and 6.69-45.51 m 3 /s for the positive response taxa (Stylonychia and Monoraphidium).

Evolution of Community Structure and Function of Bacteria and Eukaryotes
It is challenging to assess the temporal and spatial variability in complex biomes, as individual taxa respond differently to changing environmental conditions (Meng et al., 2020;Xiao et al., 2021). In spring and autumn, the water temperature varies considerably, and the types of organisms differ, resulting in different species composition of biological communities. In summer, high temperature is suitable for the growth and reproduction of thermophiles. In winter, owing to the low water temperature, biological species reduce, but suitable for the growth of some psychrophilic organisms (Salem, 2021). Due to the influence of habitat differences, nutrients, and water conservancy facilities, the biogeographic distribution of biodiversity and community structure in the upper, middle and lower reaches had significant heterogeneity (Cantonati et al., 2022). The microbial community in the middle and lower reaches is affected by more deterministic selection (Newbold et al., 2020). Previous studies have shown that microorganisms in water may be related to input from industrial sources within the watershed (Sagova-Mareckova et al., 2020). Acinetobacter is a core taxon of the sewage bacterial community and was used to indicate water contamination . Changes in hydrological factors have also led to differences in the dispersal capacity of aquatic organisms in water, enabling spatial similarity or heterogeneity in biome structure (Ćmiel et al., 2020). In addition, the restriction of biome diffusion by the Baojixia Dam might further contribute to the apparent differences in the species composition between upstream, middle, and downstream communities (Rasdi et al., 2021).
Increasing studies have highlighted the importance of bacterial and eukaryotic communities in nutrient cycling and biochemical degradation (Stewart et al., 2018). The intensity of human disturbance impacted biological communities' composition and functional attributes and the relationship between them. Previous studies indicated that nitrogen, phosphorus, and excessive nutrients had profound effects on microbial communities' functional attributes in water (Ouyang et al., 2020). We also FIGURE 10 | The correlation analysis between eco-hydrological variables and the taxa at the phylum level (*p < 0.05). Flow velocity (V), river discharge (Q), and water level (D).
Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 929541 observed that most bacteria were involved in metabolic pathways. This result is probably because the microorganisms with pollutant degradation function are selected by high-intensity environmental effects, leading to the synchronous changes of microbial community composition with the alteration of functional attributes (Dini-Andreote et al., 2015). With the increase in water temperature and flow, the microbial activity from spring to autumn was promoted, probably contributing to the degradation of pollutants by microorganisms. Consistent with previous studies, hydrochemical characteristics, temperatures, and nutrients, alone or together, affect microbial-mediated nutrient cycling (Chakraborty et al., 2021). Several studies have shown that a higher abundance of pollutant-degrading bacteria was observed in industrially polluted water, and changes in land use types promoted alterations in the metabolic functions and structures of bacterial communities (Zhang et al., 2016;Zhang et al., 2021). In this study, Cyanobacteria and Gammaproteobacteria, as the primary nitrogen-fixing microorganisms, are the dominant groups in warm seasons and downstream reaches. Previous studies have shown that nifH could be used as a characterization gene of nitrogen fixation to detect the population structure and diversity of nitrogen-fixing microorganisms in the water environment (Byappanahalli et al., 2019). The amplified sequence of nosZ gene was similar to that of nitrogen-fixing microorganisms, indicating that nitrogen-fixing bacteria are closely associated with denitrifying bacteria . The decrease in temperature and the increase of DO content would limit the function of reducing microbial bacteria in biological nitrogen removal, reducing the nitrogen source required for microbial growth (Ying et al., 2020). Most phosphorus cycle function genes with significant differences in abundance were inorganic phosphate solubilization, indicating that temperature and water quality changes mainly affect the microbial utilization of inorganic phosphorus such as phosphonates and phosphates (Luo et al., 2017). The biochemical process of phosphorus was closely coupled with the nitrogen cycle, which supports the highly significant positive correlation between the functional groups involved in the transformation of nitrogen and phosphorus (p < 0.001) (Stoliker et al., 2016).

Network Interaction Mechanism Among Multi-Trophic Groups
Most studies on river time series and geographic distribution have focused on identifying changes in limited organisms but rarely FIGURE 11 | TITAN sum (z−) and sum (z+) values corresponding to all candidate change points along the analyzed environmental gradient. Dashed and solid lines represent the cumulative frequency distribution of change points for 500 bootstrap replicates for sum z+ and sum z−, respectively. Peaks in the sums of z− and z + are the sites along the gradient with synchronic decrease and increase of taxa, respectively.
Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 929541 explored biological interactions between them (Zhao et al., 2013). The eDNA metabarcoding combined with network analysis provided a unique perspective on the interpretation of river aquatic biodiversity, offering vital information about the complex biological interactions associated with ecosystem changes by simultaneously examining multiple biological domains and trophic levels (Vasselon et al., 2017). In the molecular ecological network analysis, the modular characteristics of biological interaction networks between bacteria and eukaryotes were relatively prominent, indicating that species within the module may be more closely related in survival or function (Chun et al., 2020). The close association and modular phenomenon of taxa in the community also indicate that these biological groups occupy a similar niche and have similar ecological functions (Shi et al., 2019). Vital modules maintain connectivity in the network, conducive to defining critical species in the system. The complex interactions between multi-trophic groups may fundamentally impact species richness, community composition, ecosystem function, and stability . Multi-trophic communities are characterized by high taxonomic and functional diversity. The increase of network nodes and connectivity may provide more functional redundancy, complementary effects, and diversification effects, which was conducive to the stability of the community, so as to improve the utilization efficiency of resources and enhance the resistance and adaptability to environmental changes (Saleem et al., 2019). Heterotrophic FIGURE 12 | Taxon-specific TITAN responded to environmental change, showing significant indicator taxa (p ≤ 0.05). z + and z − represent taxa increased or decreased in the frequency of occurrence and abundance with an increased environment variable, respectively. Circles represent the change point of each taxon in proportion to the magnitude of the response.
Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 929541 bacteria, algae, protozoa (ciliophora and heterotrophic flagellates), and metazoa (arthropods, nematodes, and rotifers) constitute the main components of the microbial food web (Houliez et al., 2021). The results of this study suggested that the interaction between environmental variables and microbial network controls the relationship between different biological groups and directly affects microbial food web structure. Heterotrophic bacteria and invertebrates play a vital regulatory role in the microbiota in the microbial food web. Algae decomposition in rivers is achieved mainly through the direct feeding of invertebrates and the decomposition of extracellular enzymes derived from fungi and bacteria (Marks, 2019). Previous studies have shown that invertebrate and protozoa communities may profoundly impact the abundance, structure, and composition of primary producers, bacteria, and fungi (Yang et al., 2022). Seasonal changes affect the interaction between organisms at different nutrient levels. Due to the low temperature and the decrease of flow and velocity in winter, the enrichment of nutrients in the water body could reduce the competition of organisms for resources and the degree of niche sharing among biological groups (Saleem et al., 2019). The predator-prey relationships are essential interactions between organisms with different nutritional levels (Kruk et al., 2021). Consistent with most food web studies, we emphasized the crucial role of nutrient interactions in biota. The coexistence of bacteria, algae, and plankton in the Weihe River seasonal networks is speculated to be an interspecies predator-prey relationship. Sufficient food resources can maximize the growth rate of the microbial community and enable organisms to convert resources quickly (Brose et al., 2019). In the hubs of all ecological networks of this study, keystone species (hub) were located at the edge of the network module in different seasons. Not all OTUs in the ecological network play a dominant role in the food web and show a predator-prey relationship. Therefore, we further analyzed the relationships between microbial food web biota and its response to nutrient changes through a structural equation model. The results showed that the interaction of dominant and rare groups was vital for multi-trophic groups to maintain community relationships and enhance community stability under nutrient changes. Nutrient concentrations such as nitrogen and phosphorus are not the only variables affecting the microbial food web model. Nutrients transferred from microorganisms to animals is a vital mechanism of nutrient regeneration, especially for nitrogen and phosphorus .

Response Mechanism of Multi-Trophic Community to Eco-hydrological Changes
River ecology and hydrology changes may affect the food web model, thus altering the structure and relationship of different groups in the microbial food web (Mor et al., 2018). Previous studies have pointed out that hydrological factors play a critical role in controlling the interaction between predators and prey (Yang et al., 2015). This research reveals that the land use of farmland and urban construction reduce biological diversity and affect the community interdependence of multi-trophic groups, inhibiting typical ecosystem processes within biological communities, and reducing the ecosystem function. In addition to the direct impact of eco-hydrological variables, the indirect cascade effects may also be the reason for the notable enhancement of species interactions in microbial networks in the Weihe River. For instance, invertebrates directly affected the shaping of algae and microorganisms, thus indirectly modulating the function of ecosystems. Compared with the influence of eco-hydrological variables, top-down forces (predation/competition) also significantly impact the nutrient transfer from bacteria to invertebrates. Some studies have suggested that the process of species selection and diffusion result in the differences in the composition of heterotrophic bacteria, algae, fungi, and protozoa communities (Epsztein et al., 2020). In addition, ecological drift plays an essential role in shaping the structure of the water microbial community in the Weihe River (Huber et al., 2020).
Influenced by the hydrological regime and water environment, quantitative analysis of the response of aquatic organisms to changes in ecological and hydrological elements is crucial to river health and determining ecological flow (Xuan et al., 2019). Some studies have shown that changes in hydrology may alter plankton spatial distribution and species structures, and that temperature plays a direct role in plankton growth (Qu et al., 2018). The negative response taxa in microbial taxa, Flavobacterium, indicates that aerobic microorganisms are sensitive to environmental changes (Kim et al., 2021). In this study, there were certain but slight differences in eco-hydrological thresholds between bacteria and eukaryotes, indicating that multi-trophic groups can better determine the range of ecological thresholds (Moutinho et al., 2021). In addition, biological taxa gradually adapt and continue to survive under environmental stress (Bruno et al., 2019). The driving relationship between organisms and habitat may be responsible for the difference in indicator taxa and ecohydrological thresholds in the Weihe River (Baker and King, 2010).

CONCLUSION
In this study, eDNA technology was used to analyze the characteristics of the changing structures and functions of bacterial and eukaryotic communities. The community structures of bacteria and eukaryotes showed notable seasonal changes (p < 0.05). The multi-trophic groups in the middle and lower reaches of the Weihe River were mainly affected by habitat differences. The abundance of functional genes varied slightly (p < 0.001), among which nitrogen and phosphorus cycle function genes showed the characteristics of abundance fluctuation in spring, summer, and autumn. The interaction between biological taxa was complex, and network modularity existed within the community. Multi-trophic groups showed an apparent predator-prey relationship. In the aquatic ecosystem, a microbial food web with multi-nutrient levels played a vital role in the circulation of nitrogen and phosphorus.
According to correlation analyses, a significant association existed between bacterial and eukaryotic groups and ecohydrological variables (p < 0.05). TITAN analysis was used to identify eco-hydrological thresholds and potential indicator taxa for the multi-trophic community. Our approach demonstrated the value of multi-trophic taxa for setting goals for biodiversity conservation. However, due to the limitations of study subjects and spatiotemporal scales, we suggest that further work is required to determine suitability thresholds of taxa. The study demonstrated the application of eDNA technology in the field of river biodiversity conservation and ecological restoration. Meanwhile, it also provides a new perspective on how the food web respond to interferes with the environment and promote the nutrient cycle under the changing environment, reinforcing the concept of the assembly of bacteria and eukaryotes communities (Moutinho et al., 2021;Cantonati et al., 2022).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.