Effect of different types of industrial wastewater on the bacterial community of urban rivers

Abstract Although the health of rivers is threatened by multiple anthropogenic stressors with increasing frequency, the impacts of different types of industrial pollution with regard to the structure and diversity of the bacterial community in the river have not been fully recorded. In order to study the composition and interaction of bacterial communities and co-occurrence patterns of bacterial communities in rivers polluted by different industries, we sampled the sewage outfalls of textile factory, food manufacturer, steel plant and lighting factory. The result showed that the bacterial community composition significantly differed in the rivers polluted by distinct types of industrial pollution. In addition, correlation-based network analysis showed that the microbial under different industrial pollution types have a nonrandom modular structure, and the nodes in the different modules perform different functions. In summary, adaptive changes in the composition of bacterial communities in rivers and the occurrence of species interactions are responses to different types of industrial pollution.


Introduction
Rivers are unique ecosystems, an important part of the aquatic environment and the main water source for human consumption (Breton-Deval et al. 2019), but rapid industrial development leads to the gradual increase of pollutants in rivers (Rui and Mingming 2011), and the degradation of river aquatic ecosystems (Sheng et al. 2013;Liang et al. 2018). At present, the discharge of industrial sewage has become one of the main reasons for the deterioration of river ecology (Raat 2001;Rajaram and Das 2008), and different types of industrial pollution may break the balance of river ecological functions.
Previous studies have extensively used biota, including invertebrates and diatoms, as biological indicators of local pollution levels to evaluate the health of aquatic ecosystems (Barbour et al. 1989;Reynoldson et al. 2001;Bonada et al. 2006;Wu et al. 2013).
However, owing to their partial extinction, their diversity and abundance have been reduced, and their application in highly damaged aquatic environments has been hindered (Li et al. 2017a). Therefore, it is necessary to develop new sensitive and powerful indexes to recognize and monitor changes in the response of river ecosystems to industrial pollution.
Recently, with the expansion of high-throughput sequencing technology, people have become increasingly aware that microorganisms play a significant role in energy conversion and information transfer in aquatic ecosystems (Okafor 2011;Fan et al. 2019;Li et al. 2020) and are extremely susceptible to fluctuations in the surrounding environment (Martins et al. 2011;Oikonomou et al. 2015). Monitoring microbial community dynamics under certain ambient pressures can facilitate the identification of bacterial indicators and shed light on the evaluation and prediction of water degradation (Fortunato et al. 2013). As a microbial community in river ecosystems, the bacterial community plays an important role in nutrient cycling, energy flow and pollutant degradation (Ruggiero et al. 2006;Azam and Malfatti 2007;Inamdar et al. 2012), and the diversity of the bacterial community can reflect its ecological function and the function of river ecosystems. Thus, bacterial communities in rivers are very promising candidates for biological indicators. To date, a few studies have been carried out to provide a large amount of information on the impacts of different sources of pollution on river bacterial communities (Ibekwe et al. 2016). For instance, bacterial communities in rivers are affected by wastewater discharge from sewage treatment plants (Guo et al. 2012) and agricultural diffusion pollution (Pernet-Coudrier et al. 2012). In contrast, there are few studies on the impacts of industrial wastewater discharge on bacterial communities in rivers. Studies on the effects of different types of industrial pollution on bacterial communities in rivers are even scarcer.
Variation in bacterial community structure and diversity can be induced by various kinds of physical and chemical properties of water, including temperature, pH (Xiong et al. 2012;M endez-Garc ıa et al. 2014), nutrients (nitrogen, phosphorus) (Chrost et al. 2009;Yuan et al. 2016), metal concentrations (such as Cu, Cr, Pb) (Feris et al. 2009), and other pollutants (such as polycyclic aromatic hydrocarbons) (Sabater et al. 2016;Xie et al. 2016). These variations are frequently used to identify indicators or different bacterial species related to certain pollutants or ambient conditions. Although previous research examined changes in the composition and diversity of bacterial communities and identified bacterial indicators related to their respective environments (Wang et al. 2012;Bier et al. 2015;Kaestli et al. 2017;Richa et al. 2017), non-biological factors alone often do not fully allow the prediction of community dynamics, which means that biological factors, such as interspecific interactions, need to be incorporated (Louca et al. 2016). The interactions among microorganisms will affect the performance and function of the bacterial community and even the function of the entire ecosystem (Deng et al. 2016;Ma et al. 2016b;Hu et al. 2017;Morri€ en et al. 2017;Li et al. 2017b). Therefore, information on microbemicrobe interactions can provide new insights into the variations in bacterial communities under different environmental conditions and is helpful for identifying the corresponding bacterial indicators.
Bacterial cooccurrence networks have been used to study the interactions between microorganisms and key species in ecosystems and to explain the potential intraspecific or interspecific interactions occurring in water environments (Berry and Widder 2014). Recent studies have indicated that bacterial communities commonly have modular structures and nonrandom co-occurrence patterns, which intensely hints at the vital role of biological interactions in governing community aggregation and ecosystem function (Kara et al. 2013;Berry and Widder 2014;Cram et al. 2015). Moreover, it has been proposed that in an ecological network, some keystone species play a non-substitutable role in preserving the structure and function of the entire bacterial community (Sazima et al. 2010;Ma et al. 2016a). However, little is known about how industrial pollution affects the population dynamics of central species. Answering this significant question will enhance our understanding of the relative roles of biotic, abiotic and anthropogenic factors in bacterial composition, diversity and species interactions and will shed light on the ecological rules that guide bacterial community assembly in river ecosystems. Simultaneously, functional prediction technology can provide a glimpse of the general picture of the functional spectrum of bacteria, fully utilizing the advantage of the highly cost-effective sequencing of the diversity composition spectrum of bacteria to realize the prediction of the metabolic function of bacteria (Mehetre et al. 2016). Here, we used high-throughput Illumina MiSeq sequencing technology to analyze the characteristics of a bacterial community through the combination of microbe-environment correlations and microbe-microbe interactions to determine whether under four different types of industrial pollution (light industry, manufacturing industry, modern industry, and heavy industry) (1) there are differences in the bacterial communities in rivers; (2) the co-occurrence pattern of river bacterial communities is random and driven by certain factors; (3) there are indicator species in river bacterial communities that can maintain the structure and function of the ecological community.

Sites description
The Qingliu River is a typical urban river in Chuzhou City, Anhui Province, China (with a drainage area of 1252 km 2 and a total length of 84 km) and is the left-bank tributary of the Chuhe River, a first-order tributary of the Yangtze River. With the rapid increase in population density, river pollution has become an intensifying problem. The large-scale discharge of different types of industrial wastewater along the river is one of the important causes of aquatic environment pollution. In this study, four unique sampling points were selected from the upper reaches to the lower reaches of the Qingliu River, which were located at the sewage outlets of a textile mill (FZW), steel plant (GGW), food factory (SPW) and lighting factory (ZMW) (Figure 1). These four sampling points cover manufacturing, heavy industry, light industry and modern industry.

Sample collection and physicochemical analysis
Sampling was carried out in May 2019. Water samples (approximately 0.5 m sampling depth) were collected with a plexiglass water sampler, and the samples from each river were collected in triplicate. The three subsamples from each site were randomly obtained at 3 m intervals. The pH, temperature (T) and dissolved oxygen (DO) of the water were determined immediately after sampling with a multiparameter controller (YSI 6600V2, USA). All the samples were shipped within 1 hour to the laboratory in a cooler and then kept at 4 C for further analysis within 2 days.
The NH 4 þ -N (ammonium), TP (total phosphorus) and TN (total nitrogen) in the water were measured in the laboratory according to standard methods (Jin and Tu 1990). The content of metal elements Cu, Zn, Pb, Cr, Cd in water samples were determined by using the X7 inductively coupled plasma mass spectrometer (ICP-MS) of Thermo Corporation in the collision cell mode.

DNA extraction
Using a FastDNA Spin Kit (MP Biomedicals, USA) to extract DNA from each water sample within one day after sampling on the basis of the manufacturer's instructions. The water samples were filtered through a 0.2 lm filter and finally cut into small pieces for DNA extraction. DNA quality was checked on a 1.5% agarose gel. Using an ultraviolet spectrophotometer (Eppendorf, Germany) to quantify DNA. DNA was stored at À20 C until further processing (Keshri et al. 2018).

16s rRNA gene analysis and sequence analysis
The universal primers 338 F (5 0 -ACTCCTACGGGAGGCAGCAG-3 0 ') and 806 R (5 0 -GGACTACHVGGGTWTCTAAT-3 0 ) were used to amplify the V3 to V4 regions of the bacterial 16S rRNA gene. The PCR condition was as follows: 94 C for 5 min, 30 cycles of 94 C for 30 s, 54 C for 30 s, 72 C for 45 s, and finally 72 C for 10 minutes (Chen et al. 2017). Process all samples on the basis of the formal experimental conditions and repeat each sample 3 times. An AxyPrep DNA Gel Recovery Kit (AXYGEN Company) was used to cut the gel to recover the PCR products. The Tris-HCl buffer was used for elution and 2% agarose was detected by electrophoresis. With reference to the preliminary quantitative results of the electrophoresis, the PCR products were analyzed by a QuantiFluorTM-ST Blue Fluorescence Quantitative System (Promega Company) perform detection and quantification, and then mix the corresponding proportion according to the sequencing quantity requirements of each sample. Finally, sodium hydroxide was used to denature to produce single-stranded DNA fragments. Using Illumina Miseq system to perform 16S ribosomal RNA gene sequencing (Illumina Inc., San Diego, CA, USA). Using Trimmomatic to demultiplex and filter raw pyrosequencing data (Bolger et al. 2014). FLASH (v1.2.11) (Mago c and Salzberg 2011) was used to amalgamate the overlapping reads into a single long read (Lie et al. 2014). VCHIARCH (v1.11.1) (Rognes et al. 2016), and UCHIME (v7.0 http://drive5.com/uparse/ ) algorithm (Edgar et al. 2011) were used to examine potential chimeric sequences. According to different similarity, all sequences are divided to clustered, and a statistical analysis of 97% similarity is usually performed on OTUs (Schmidt et al. 2013). Use the RDP classifier for obtaining the species classification information corresponding to each OTU (v2.11 https://sourceforge.net/projects/rdp-classifier/). Bayesian algorithm (confidence threshold 0.7) on 97% of OTU Classification was performed using the QIIME (v1.9.1) platform (http://qiime.org/install/index.html) to represent sequences at similar levels (Samarajeewa et al. 2015), and the community composition of each sample was analyzed. The raw sequencing data have been submitted to NCBI Sequence Read Archive under the project accession code PRJNA600788.

Network analysis and statistical analyses
Using the R package Hmisc pair to calculate the Spearman rank coefficient (p) among all 97%-cutoff OTUs at least 12 OTUs occurrence and at least 10 abundances in each group (Harrell and Dupont 2008, http://cranr-projectorg/web/packages/Hmisc/indexhtml). Pick out OTUs with robust (p ! 0.8) and significant (P value < 0.01) correlations. Subsequently, network visualization and modular analyses were carried out on using Gephi (http://gephi. github.io/). Gephi determined the average path length (the average number of steps along the shortest path for all probable network node pairs) and the topological properties of each node including degree (number of adjacent edges) and betweenness centrality (number of shortest paths through nodes). Defined nodes with low betweenness centrality values and high degree as keystone species in co-occurrence networks (Berry and Widder 2014), which is similar to the central species used in other studies (Murphy et al. 2014). Each node size is directly proportional to the relative abundance of the bacterial community (1og (n þ 1)).
Statistical analysis was performed using R3.5.1, including mean value, standard deviation, Kruskal-Wallis rank tests (KW) and Spearman correlations. It is reported that t-test results with p 0.05 were significant. Using redundancy analysis (RDA) in R software (v3.5.1, vegan package) evaluated correlations between environmental variables and bacterial community structure.

Bacterial community diversity
In this study, we used Illumina MiSeq high-throughput sequencing to analyze the bacterial community diversity, composition and cooccurrence patterns in twelve samples of wastewater from a textile mill, steel mill, food plant and lighting plant. After the removal of low-quality sequences and chimeras, a total of 359,028 effective 16S rRNA sequences and 601 OTUs were acquired. All the sparse curves tend to approach the saturation plateau, and the sequencing data reach saturation, indicating that the sequencing depth is sufficient for community analysis. The coverage of all samples averaged 99.79 ± 0.05%. The good coverage value of the samples shows that the distribution of the bacterial flora in the four industrial wastewaters is correct. In addition, the number of OTUs in the 12 samples from four factories differed. Among them, the OTUs of FZW and GGW were 356 and 276, respectively. The OTUs of SPW and ZMW were 463 and 443, respectively, and were similar. The results indicate that SPW has the highest species richness and that GGW has the lowest species richness.
The Shannon index reflects bacterial diversity. FZW and GGW had values of 3.71 and 2.82, respectively, while SPW and ZMW had values of 4.18 and 4.08, respectively. Hence, the diversity of SPW and ZMW is higher than that of FZW and GGW. The Simpson index reflects the uniformity of bacterial diversity, with the largest value being observed for GGW, the smallest for SPW, and FZW and ZMW having values of 0.0879 and 0.0726, respectively ( Figure 2; Table 1). The results indicate that GGW has the highest evenness component of diversity.

Bacterial community structure
The beta diversity of the 12 industrial wastewater samples was studied according to PCoA based on the Bray-Curtis distance ( Figure S1(a)). Generally, the closer two samples are plotted, the more similar their compositions are. It is worth noting that SPW and ZMW had similar microbial communities, while the samples from FZW were marginally different from the others. However, the bacterial composition of GGW was different from that of the other three samples. ANOSIM (r ¼ 0.89, P ¼ 0.001) indicated that the four sampling points had significantly different bacterial communities ( Figure S1(b)).
Further analysis indicated that class Alphaproteobacteria, class Betaproteobacteria, and class Gammaproteobacteria, members of Proteobacteria, exist in the samples. It is worth noting that class Gammaproteobacteria was dominate in GGW (Figure 4(a), Kruskal-Wallis tests, p < 0.05). Relative to that at the other sampling points, the relative  abundance of the family Mycobacterium in FZW was obviously higher at the family level (Figure 4(b), Kruskal-Wallis tests, p < 0.05). In addition, the classification OTUs and relative abundances were analyzed at the genus level (Figure 3(b)). Some genera appeared in each sample but were more abundant in individual samples than in others, such as the genus Acinetobacter (Figure 3(b)). The relative abundance of the genus Acinetobacter in GGW was higher than that at the other sampling points (45.10 ± 38.77%). Genus Exiguobacterium had a relative abundance (8.12 ± 5.72%) at ZMW that was higher than that at the other sampling points. The genus Flavobacteria was significantly higher at SPW than at the other sampling points (Figure 3(b), Kruskal-Wallis tests, P < 0.05). The unique genus at GGW is the genus Bifidobacterium (1.10 ± 0.95%). The difference in the species and quantity of bacterial genera in the industrial wastewater in the different treatment units indicates that the distribution of bacterial genera in the samples is diverse.

Multivariate correlation analysis between water bacterial communities and environmental variables
Environmental factor analysis compares microorganisms with environmental variables according to their classification level and evaluates the correlation between microorganisms and environmental variables. The environmental parameters at the different sampling points are shown in Table S1, and the physical and chemical indexes of the 12 samples are different. Except at FZW, the pH value of the water samples was between 7.24 and 7.53, with a small range of variation, while the pH value at FZW was 9.96, which is slightly alkaline. The temperature of the water samples showed little variation, ranging from 21.2 C to 22.7 C. The TP, TN, NH 4 þ -N and DO at SPW were higher than those at the other sampling points (Table S1). Furthermore, to understand the correlation between the composition of the microbial communities and environmental factors in rivers polluted by different industries, RDA was performed ( Figure 5(a)). Among all the environmental factors examined, pH (r 2 ¼ 0.730, P ¼ 0.003) was significantly related to the microbial structure in the water, while the correlations between T (r 2 ¼ 0.075, p ¼ 0.806), TP (r 2 ¼ 0.109, p ¼ 0.611), TN (r 2 ¼ 0.417, p ¼ 0.122), NH 4 þ -N (r 2 ¼ 0.215, p ¼ 0.355) and DO (r 2 ¼ 0.136, p ¼ 0.546) and microbial structure was not significant.
The results show that pH is the most important environmental factor affecting the water samples. In addition, in the water samples, 37.65% of the constrained variance was explained by RDA1, and 35.66% was explained by RDA2 ( Figure 5(a)).

The topological and taxonomic properties of the co-occurrence networks
In light of the nonrandom community assembly model in the river bacterial community, a network interface was constructed to further probe the topological and classification characteristics of the bacterial co-occurrence model (Figure 6(a)). In view of the correlation analysis, the 4,015 edges represent a very strong cooccurrence correlation among the 189 OTU nodes (r ! 0.8, P < 0.05), and the characteristics of the network are that the positive correlation is much higher than the negative correlation. The networks were represented by a modular value (MD ¼ 0.658), average path length (APL ¼ 2.154) and average degree (AD ¼ 42.487) that were higher than those of other global networks, indicating that the distinct types of industrial pollution resulted in apparent cooccurrence patterns and modular structures. Modularity analysis indicated that the whole network can fall into 6 main modules. Each main module consists of a group of OTU nodes. Compared with that of the nodes in other modules, the interconnection frequency between them is higher (Figure 6(b)). OTUs from module I had the highest relative abundances in GGW. Module II and module IV had the highest relative abundances in SPW, while module III and module V commonly occurred in ZMW and FZW. In addition, taxonomic relatedness is obviously the critical factor in determining the modular structure of the network (Figure 6(b)). For instance, the genus Acinetobacter and the genus Brevundimonas have advantages in module I (10%) and module III (8.6%) respectively. The genus Flavobacterium was the dominant genus (13.3% and 17.8%) in module II and module IV, respectively, while the genus Mycobacterium was the dominant genus (11.8%) in module V.
Keystone species play a nonsubstitutable role in maintaining the structure and function of a bacterial community. In the co-occurrence networks, they were considered to have a high degree (>70) and low betweenness centrality values (<60) (Ma et al. 2016a). According to this standard, 16 OTUs were recognized as keystone species, accounting for 8.5% of the total reads. The taxonomic classification indicated that the major keystone species belonged to the genus Candidatus Rhodoluna (3 OTUs and 14.3% in the central community), genus Pseudarcicella (1 OTUs/7.1%), genus Candidatus Planktophila (1 OTU/7.1%) and genus Flavobacterium (1 OTU/7.1%).

Bacterial metabolism
Bacterial functional genes involved in carbon (C) cycling According to the KEGG enzyme profiles, the major functional genes related to C cycling were identified from all water samples in this research. Gene abundances in active C degradation in water samples include pullulanase and alpha-amylase for starch decomposition, endo-1,4-beta-xylanase and alpha-N-arabinofuranosidase for hemicellulose decomposition, and 2-isopropylmalate synthase, limonene-1,2-epoxide hydrolase (dimEH) and isocitrate lyase for decomposition of aromatic substances. In all samples, some functional genes, such as RuBisCO (ribulose-bisphosphate carboxylase small chain), glyoxal oxidase, lignin peroxidase and cellobiase related to intractable C degradation were absent. Genes related to functioning of C fixation such as acetyl-/propionyl-CoA carboxylase, carbon monoxide dehydrogenase/acetyl-CoA synthase subunit alpha, transketolase, glyceraldehyde-3-phosphate dehydrogenase (NADP), fructose-bisphosphate aldolase, pyruvate kinase isozymes R/L, 2-phosphoglycerate kinase, S-adenosylmethionine synthetase, Figure 6. Co-occurrence networks of 97%-cutoff OTUs within lotic microbial communities based on correlation analysis. A connection stands for a very strong (Spearman's q ! 0.8) and significant (FDR-adjusted P-value < 0.01) correlation. The size of each node is proportional to the number of connections (i.e. degree). (a) OTUs colored by phylumlevel taxonomy (b) OTUs colored by modularity.
anthranilate synthase/phosphoribosyl transferase and phosphoribulokinase were present in abundant numbers in all water samples.

Bacterial functional genes in nitrogen (N) cycling
Genes responsible for ammonification included urease subunit gamma/beta and glutamate dehydrogenase (NAD (P) þ); for assimilatory N reduction included nitrite reductase (NAD(P)H) small subunit and ferredoxin-nitrate reductase; for denitrification included nitric oxide reductase FlRd-NAD (þ) reductase and nitrous-oxide reductase; for dissimilatory N reduction to ammonium incduded periplasmic nitrate reductase; for nitrification included nitrite reductase (NAD(P)H) small subunit for nitrifiers(an indication of nitrificaion activity) and ammonia monooxygenase subunit C; and for nitrogen fixation included nitrogenase (flavodoxin).
Bacterial functional genes involved in phosphorus (P) and sulfur (S) cycling Genes for P utilization included polyphosphate kinase, exopolyphosphatase/guanosine-5 0triphosphate, 3-phytase and 4-phytase/acid phosphatase. The abundant genes for adenylyl sulfate reductase, subunit B, and sulfite reductase, dissimilatory-type beta subunit for sulfite reduction and sulfur oxidation, respectively, were detected in all water samples, ensured bacterial role in sulfur cycles.

Discussion
Different types of industrial pollution may break the balance of river ecological functions, and bacterial communities can effectively promote the circulation of organic matter (Ruggiero et al. 2006;Inamdar et al. 2012), which plays a vital role in the environment and ecological processes of river ecosystems (Araya et al. 2003). In this study, the bacterial communities in rivers affected by different industrial pollutants were studied by highthroughput Illumina MiSeq sequencing technology. The main phyla in this study were Proteobacteria (57.42%), Actinomycetes (16.75%), and Firmicutes (8.13%) (Figure 3(a)). This is consistent with findings found in other studies (Staley et al. 2013;Bai et al. 2014;Ruiz-Gonz alez et al. 2015;Ibekwe et al. 2016).
As the most dominant phylum at all sampling points, Proteobacteria is usually one of the most important prokaryotes in the water (Garrity et al. 2015), and most of its taxonomic groups are involved in various biogeochemical processes in the aquatic ecosystem Zhang et al. 2014). Further analysis showed that, at the class level, the proteobacteria in the sample were divided into Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria. In contrast to the other samples, Gammaproteobacteria were dominant in the GGW samples. Gammaproteobacteria is regarded as the most abundant group in the sediments of eutrophic lakes (Bai et al. 2012), and the extreme enrichment of these bacteria in association with iron and steel plants may be due to the organic emissions from iron and steel plants. Studies have shown that there is a remarkable negative correlation between Gammaproteobacteria and NH 4 þ -N (Zhang et al. 2015). However, in this study, the RDA indicated a remarkable negative correlation between Gammaproteobacteria and TN ( Figure 5(b)), which indicates that Gammaproteobacteria may be inseparable from nitrogen transformation and cycling in aquatic ecosystems. In addition, the genus Acinetobacter, a genus of Proteobacteria, was widely enriched in ZMW and GGW and can be used to bioassay the ecotoxicity of wastewater contaminated with heavy metals (Abd et al. 2006). Moreover, some studies have indicated that wastewater containing heavy metals from iron and steel plants can decompose soil and release heavy metals into rivers through acidic runoff (Holding 2004), which may be an important reason for the extensive enrichment of the genus Acinetobacter in association with lighting plants and iron and steel plants found in this study. There is no doubt that the abundance of this genus can be used as a sensitive indicator of the content of heavy metals in wastewater discharged from lighting plants and iron and steel plants into rivers.
As the second most enriched phylum in this study, Actinobacteria play a vital role in the carbon cycle in freshwater ecosystems (Rifaat 2003). The RDA showed that there was a remarkable positive correlation between Actinobacteria and TN and TP (Figure 5(a)), which also suggests that Actinobacteria are active in eutrophic environments. Compared with that at the other sampling sites, the relative abundance of the family Mycobacterium in Actinobacteria was higher at FZW. Family Mycobacterium (Actinobacteria) is able to mineralize some refractory polycyclic aromatic hydrocarbon (PAH) compounds with low water solubility and low bioavailability (Jazestani and Prasher 2011). However, in the textile industry, polycyclic aromatic hydrocarbons (PAHs) are mainly discharged into rivers as the basic raw material of synthetic dyes. At present, many industrial sites have been contaminated with toxic hydrophobic organic pollutants (HOCs), and PAH is one of the culprits (Liu et al. 2007). Thus, detecting the relative abundance of the family Mycobacterium in a river can effectively indicate the level of HOCs in the river. In addition, further studies have found that the genus Bifidobacterium (Actinobacteria) is a unique genus at GGW (Figure 3(b)). Undoubtedly, there are many complex metal ions in the wastewater from iron and steel plants, such as Fe 2þ , Cu 2þ , Pb 2þ and Cd 2þ (Table S2), and the genus Bifidobacterium (Actinobacteria) can take up iron in the presence of Ca 2þ and Mg 2þ Bezkorovainy 1991, 1993), indicating that the genus Bifidobacterium plays an irreplaceable role in the detection of iron in the wastewater from iron and steel plants.
It is worth mentioning that the abundance of Bacteroides at SPW is significantly higher than that at the other sampling points (Figure 3(a)). Bacteroidetes is a member of Bacteroides with strong phenotypic and metabolic diversity (Tanaka et al. 2017), and has been identified as a denitrifying agen (Liu et al. 2008). In addition, members of the genus Flavobacteria, belonging to Bacteroidetes, have a relatively high abundance at SPW (Figure 3(b)) and are especially common and abundant heterotrophs in fresh water (Michaud et al. 2012;Garc ıa-Armisen et al. 2014), playing a key role in the biogeochemical processes of the cycling of carbon and other elements (Kirchman et al. 2004;Teeling et al. 2012;Williams et al. 2013). Therefore, we believe that this genus is an important indicator for the detection of organic substances in food plant wastewater. Our results indicate that the bacterial communities in rivers are driven by various environmental factors, and different types of industrial pollution have significantly different effects on the diversity and composition of bacterial communities in the same river.
The topological characteristics of the network normally reflect the interactions among microorganisms. The betweenness centrality indicates the potential of control exerted by a single node on the interactions of other nodes in the network and emphasizes its usefulness in defining the key species in the system (Vick-Majors et al. 2014;Banerjee et al. 2016), while the degree value is a local quantification characteristic of the direct co-occurrence interactions of particular nodes (Whitman 2015). In this study, we identified the genus Candidatus Rhodoluna and genus Pseudarcicella as the two keystone taxa according to their high degree (> 70) and low betweenness centrality values (< 60). The genus Candidatus Rhodoluna is affiliated with the family Microbacteriaceae of the phylum Actinobacteria, and its members are able to assimilate ammonia (Greenblum et al. 2012), and the genus Pseudarcicella is implicated in lowering the nitrate content, which can be reduced to nitrite (Horsley 1979), and is mainly found at SPW (Figure 3(b)). This indicates that the nitrate and ammonia nitrogen contained in the wastewater from the food factory have obvious driving effects on the genus Candidatus Rhodoluna and the genus Pseudarcicella. Furthermore, these two main taxa play a vital role in defending the structure and function of ecological communities.
By integrating the OTU nodes of the taxonomic assignments and network structure, we found that the assembly of the bacterial communities was determined by taxonomic correlations and nonrandom, i.e. intimately related taxa were inclined to be interconnected and clustered together (Figure 6(a)) (Hu et al. 2017). In this study, our network nodes are divided into six major modules, and previous studies have shown that nodes in different modules perform different functions (Newman 2006). In the network of module I and module III, some bacteria are related to the detection, absorption and resistance of heavy metals. For example, in module I, the genus Acinetobacter (Figure 3(b)), which is dominant in GGW, can be used to monitor the toxicity of heavy metals in wastewater (Abd et al. 2006), while the genus Bifidobacterium can ingest iron under certain conditions (Kot and Bezkorovainy 1993). In module III, the dominant genera at ZMW, genus Brevundimonas and genus Chryseobacterium, have certain resistance to the bioaccumulation capacity of the heavy metals Pb, Hg, and Zn (Benmalek et al. 2014;Hamzah et al. 2015). Their dominance in module I and module III showed that heavy metal ions (Table  S2) in wastewater from iron and steel plants and lighting plants are the key factors driving these bacteria. The genus Flavobacterium is relatively abundant at SPW (Figure 3(b)), and it is the dominant genus in module II and module IV. The genus Flavobacterium can reduce nitrate and plays an important role in denitrification (Horsley 1979). There is no doubt that nitrate and other substances in the wastewater of food factories have a significant impact on the function of the genus Flavobacterium in modular structure. The genus Mycobacterium and genus Acidovorax, which are significantly enriched at FZW and dominate in module V (Figure 3(a)), contain many humus environmental isolates that can degrade many organic pollutants, including polycyclic aromatic hydrocarbons (PAHs) (Kim et al. 2010;Singleton et al. 2013). It is shown that the organic pollutants present in the wastewater from the textile mill are multifunctional compounds that manipulate the roles of the genus Mycobacterium and genus Acidovorax in the modular structure. Therefore, the complexity of community structure and functional processes following pollution by different types of industry in the same river and the co-occurrence mode of the bacterial community are nonrandom and function driven, which shows that the role of deterministic processes in modular structure is obvious, which is consistent with previous research (Jiao et al. 2016).

Conclusions
In this study, 16S rRNA gene sequencing technology and network analysis was used to study the response mechanism of bacterial community composition, diversity and cooccurrence pattern in river water to different industrial pollution inputs. Our results showed that the bacterial community in the river is driven by various environmental factors, and different types of industrial pollution had different effects on the composition and diversity of the bacterial community in the same river. In addition, the complexity of the community structure and functional processes following contamination by different types of industry in the same river and the co-occurrence pattern of the bacterial community are nonrandom and function driven, which shows that the role of deterministic process in modular structure is obvious. Meanwhile, a large number of functional genes related to C, N and P cycles were found in various river samples. The results of research on the diversity of bacterial communities in rivers under the influence of different types of industrial pollution provide a basis for monitoring and controlling the environmental pollution of water. At the same time, it provides representative data for the comparison between urban rivers and other rivers, which helps people to understand the microbial ecology of urban fresh water.

Disclosure statement
The authors declare no conflict of interest.