Interactome‐based abiotic and biotic impacts on biodiversity of plankton communities in disturbed wetlands

Despite great efforts for conservation, biodiversity in wetland ecosystems is still losing at an alarming rate. Thus, it is crucial to deeply understand ecological processes and mechanisms that potentially affect the loss of biodiversity.


| INTRODUC TI ON
Benefiting from the high productivity and distinct hydrological features, wetlands are considered as one of the most productive ecosystems on the Earth (Gibbs, 2000;Tockner & Stanford, 2002).
With less than 3% of the Earth's surface, wetlands provide habitats for ~20% of known species and contribute more than 40% of the ecosystem services globally (Fang et al., 2006;Finlayson, Davidson, Spiers, & Stevenson, 1999;Zedler & Kercher, 2005). It has been widely acknowledged that the high biodiversity is essential to support and maintain ecosystem services and functioning of wetlands (Fang et al., 2006;Meli, Benayas, Balvanera, & Ramos, 2014;Myers, 1996). However, biodiversity in wetlands, as well as wetland itself, is losing at an unprecedented rate, mainly due to a wide variety of anthropogenic disturbances Erwin, 2009;Hu, Niu, Chen, Li, & Zhang, 2017;Zedler, 2000). Since 1900, 64%-71% of the world's wetlands have been lost, with a much faster rate during the 20th and early 21st centuries (Davidson, 2014). As a result, severe loss of biodiversity has been reported along with the rapid degradation and loss of wetlands (Davidson, 2014;Fang et al., 2006;Gibbs, 2000). Therefore, effective strategies are urgently needed to conserve and restore wetlands, particularly biodiversity in wetlands.
One crucial prerequisite for conservation is to understand crucial ecological processes, as well as influential factors, that drive the degradation and loss of wetland biodiversity (Chaparro, Horváth, O'Farrell, Ptacnik, & Hein, 2018;Xu et al., 2019). However, most studies focused on the identification of environmental stressors underlying the degradation and loss of wetlands, while paid limited attention to factors that were associated with biodiversity (Meli et al., 2014;Meng et al., 2017). Although threats faced by wetlands could (in)directly affect biodiversity, the investigation of major ecological processes and associated mechanisms has been highly obstructed by complex species-environment and species-species interactions (Bortolotti, Vinebrooke, & Louis, 2016;Mensing, Galatowitsch, & Tester, 1998;Soininen, Kokocinski, Estlander, Kotanen, & Heino, 2007). Such a challenge is further complicated by the increasing anthropogenic disturbances, which could lead to severe fragmentation and excessive heterogeneity of habitats at fine geographical scales, even before significant patterns of large-scale wetland transformation are observed (Fang et al., 2006;Gibbs, 2000;Wang, Fraser, & Chen, 2017;Yang, Deng, & Cao, 2016). All of these challenges call for attention towards the understanding of mechanisms that affect biodiversity in wetlands. As recommended by Zedler (2000), community succession theories can effectively guide practical biodiversity restoration. Thus, comprehensive investigation of factors that influence community succession would benefit the understanding of biodiversity loss, identifying conservation priorities and improving the predictability of conservation programmes.
A majority of literature have documented the influence of environmental factors on wetland biodiversity, such as nutrient enrichment derived from human activities (Fang et al., 2006;Meng et al., 2017;Tockner, Pusch, Borchardt, & Lorang, 2010). However, recent studies suggest that abiotic impacts alone seem insufficient to explain the biodiversity variation, emphasizing the necessity to investigate biotic impacts (e.g., inter-and intraspecific interactions; Lima-Mendez et al., 2015). Indeed, biodiversity is a combination of biological entities and processes, rather than a simple collection of individual species (Jordano, 2016;Wang & Brose, 2018). Species in a community are involved by a myriad of complex interactions, in the form of predation, competition, mutualism, symbiosis, parasitism and others (Jordano, 2016;Zhou et al., 2010). Thus, species interactions form the architecture of biodiversity and determine ecosystem functioning (Bartomeus et al., 2016). Regarding the high biodiversity in wetland ecosystems, it is reasonable to speculate the existence of highly complex inter-and intraspecific interactions. However, previous studies of wetland biodiversity mainly focused on one single group or several separated groups, while rarely considering biotic interactions (Zhang et al., 2018). One important reason is the impediment to characterize and analyse the complex interactions. Fortunately, with the development of multi-omic techniques and bioinformatic methods, the interactome-based methods have allowed to discover complex species interactions (Bartomeus et al., 2016;Faust, Lahti, Gonze, De Vos, & Raes, 2015). Based on the interactome concept, the complexity of interactions could be well represented and modelled as networks (Layeghifard, Hwang, & Guttman, 2017). A key feature of network theory is that the architectural features of networks appear to be universal to most complex systems (i.e., scale-free), thus paving the way to develop and characterize networks involving both biological and non-biological components (Layeghifard et al., 2017). Henceforth, interactome-based methods have proven to be useful in resolving the succession dynamics of microbiome systems or the marine plankton communities (mainly zooplankton and phytoplankton; Guidi et al., 2016;Layeghifard et al., 2017;Lima-Mendez et al., 2015).
Sanjiang wetland is a promising ecosystem to investigate the biotic interactions and abiotic impacts on biodiversity in wetlands. Sanjiang wetland is located at the Sanjiang floodplain in north-east China, covering an area of 108,900 km 2 . Sanjiang Plain is formed by inundation and exposure of three large rivers-Wusuli, Songhua and Heilong Rivers, and thus providing highly diverse types of habitats. A rich biota has been recorded in Sanjiang wetland, especially endangered and vulnerable waterbirds such as Oriental Stork . Similar to many wetlands, Sanjiang wetland has been severely degraded and lost, especially in the past five decades. It has been estimated that about 70% of natural wetlands were lost from 1976 to 2005, with more than 90% of area being transformed to croplands (Song et al., 2014). The degradation and loss of wetlands have resulted in the shrinkage of suitable habitats and significant decline of biodiversity . In order to facilitate restoration and conservation programmes by the Chinese government (Lu et al., 2016), it is urgent to improve the understanding of processes and mechanisms responsible for biodiversity variability, especially under the severe pressure of anthropogenic disturbances.
In this study, we sampled two types of plankton communities (i.e., bacterioplankton and zooplankton) to investigate the composition and geographical distribution patterns of biodiversity and associated mechanisms responsible for observed patterns in Sanjiang wetland. These plankton communities were chosen as they play vital roles in wetland energy flow, material cycle and ecosystem functioning (Jansson, Andersson, Berggren, & Leonardson, 1994;Soininen et al., 2007). We sampled five locations with different levels of anthropogenic disturbances and habitat types. We characterized the biodiversity of plankton communities by high-throughput sequencing-based metabarcoding and analysed the complex biological interactions based on the interactome concept. We aimed to answer three questions: (a) Would the biodiversity of plankton communities show variability in respect to different habitats at fine geographical scales? (b) What factors drove the variability of community structure and spatial patterns of biodiversity within and among habitats? (c) What was the contribution of biotic and abiotic factors to the observed patterns?

| Collection of biological samples and measurement of physicochemical variables
Samples of bacterioplankton and zooplankton were collected from five habitats in Sanjiang wetland (Figure 1), including one marsh Wetland habitat (MW), one lacustrine Wetland habitat (LW), one paddy Wetland habitat (PW) and two riparian Wetland (RW) habitats. The LW habitat is located at the core area of Sanjiang Wetland National Nature Reserve and is better protected than the other habitats . The PW habitat has recently been restored to wetland after decades of agricultural farming. In total, 28 sites were collected, including four, five, five, seven and seven sites in MW, LW, PW, RW1 and RW2, respectively. Geographical location of each site was recorded using a Handheld Garmin GPS Navigator. A total of 500 ml water was collected at each site and filtered through a 0.22-μm membrane to collect microbe cells. Membranes were stored at −80°C until further analyses. A total of 60 L water was filtered through a 35-μm-mesh net to collect zooplankton samples, and all collected animals were immediately preserved in 100% alcohol with a final volume of 100 ml. Another 500 ml water was collected for water physicochemical analysis. Water and biological samples were stored at 4°C before further analyses.
In the field, water quality parameters were measured, including water temperature (T), electric conductivity (EC), pH, oxidationreduction potential (ORP) and total dissolved solids (TDS) using a multi-parameter water quality sonde (MYRON Company, USA).
Dissolved oxygen (DO) and concentration of chlorophyll a (Chl_a) were measured using a portable dissolved oxygen meter (Hach Company, USA) and a Handheld Fluorometer (Turner Designs, USA), respectively. The alkaline potassium persulfate digestion-based UV spectrophotometry, Nessler's reagent-based spectrophotometry and differential UV spectrophotometry were employed to measure total nitrogen (TN), ammonia (NH 4 ) and nitrate (NO 3 ), respectively.
The potassium persulfate digestion-based UV spectrophotometry and ammonium molybdate-based spectrophotometry were employed to measure total phosphorous (TP) and soluble reactive phosphorous (SRP). Water samples were pre-processed with digestion kit (HACH Company, USA) and measured using UV spectrophotometry to calculate the concentration of chemical oxygen demand (COD).
The concentration of total organic carbon (TOC) was measured using Shimadzu-TC (Shimadzu, Japan). Concentration of potassium (K), calcium (Ca), sodium (Na) and magnesium (Mg) was measured F I G U R E 1 Sampling area and sites in Sanjiang wetland. LW, lacustrine wetland (blue dots); MW, marsh wetland (green dots); PW, paddy wetland (red dots); RW, riparian wetlands (grey dots and open circles) using the inductively coupled plasma optical emission spectrometry.

| DNA extraction, PCR amplification and highthroughput sequencing
A modified CTAB extraction protocol was employed to extract total genomic DNA of bacteria (Huang et al., 2009;Yang et al., 2016). DNA extracts were then used as PCR templates to amplify the hypervariable V4 region of 16S rDNA with the primer pair of 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT; Caporaso et al., 2011). Primers were modified with an addition of 12-nucleotide tag at the 5′-end to distinguish samples. The PCR amplification program was 4 min of denaturation at 94°C, 35 cycles of 45 s at 94°C, 30 s at 50°C and 90 s at 72°C, followed by a final elongation for 10 min at 72°C. The high-throughput sequencing was performed using the Illumina Miseq PE250 platform.
Total genomic DNA of zooplankton was extracted using the DNeasy Blood and Tissue Kit (Qiagen). DNA extracts were then used as PCR templates to amplify the V4 region of 18S rDNA with the primer pair of Uni18S/Uni18SR, which was specifically designed for zooplankton communities (Zhan et al., 2013). Primer for each sample was labelled with an addition of a unique eight nucleotide tag at the 5′-end to allow pooling all samples together. The PCR amplification program was 5 min of denaturation at 95°C, 35 cycles of 30 s at 95°C, 45 s at 50°C and 90 s at 72°C, followed by a final elongation for 10 min at 72°C. The high-throughput sequencing was performed using the Illumina Miseq PE300 platform.

| Bioinformatics analyses
Raw data of bacterioplankton were sorted into independent files based on unique tags. Paired-end reads were merged before trimming tags and primers. We discarded low-quality sequences: (a) containing "N"s (undetermined nucleotides); (b) with quality score less than 20 (<Q20); and (c) with length shorter than 200 bp (Yang et al., 2016. Filtered sequences were clustered into operational taxonomic units (OTUs) at the similarity threshold of 97% using UPARSE pipeline (Edgar, 2013). Using Ribosomal Database Project (RDP) classifier, taxonomic assignment of representative sequences was obtained by searching against the SILVA_128 database at the 85% confidence level (Wang, Garrity, Tiedje, & Cole, 2007). OTUs classified as unknown, archaea and mitochondria were removed from further analyses. To improve normality of sequence data, the OTU table was rarefied down to 19,626 sequences per sample according to the smallest value of total number of sequences. All analyses were conducted using the in-house pipeline Galaxy (http://mem. rcees.ac.cn:8080/root; Feng et al., 2017).
Raw data of zooplankton were processed using the UPARSE algorithm embedded in USEARCH (Edgar, 2010(Edgar, , 2013. Artificial primers and tags were trimmed with Python scripts. We removed sequences: (a) that contained any undetermined nucleotides ("N"s); (b) with quality score lower than 20 (<Q20); and (c) that had the maximum expected error threshold lower than 0.75 (Yang et al., 2016(Yang et al., , 2018. Filtered sequences were trimmed to the same length of 225 bp based on the recommendation of Xiong et al. (2017). OTUs were clustered at a similarity threshold of 97% with de-replicated sequences Zhan, Bailey, Heath, & MacIsaac, 2014). By searching against the online nucleotide database of GenBank, taxonomic information of OTUs was obtained. All OTUs and representative sequences were filtered with parameters of e-value < 10 −80 , minimum query coverage > 80% and similarity > 85% (Zhan et al., 2014). OTUs assigned as metazoan zooplankton were kept for further analyses. For both bacterioplankton and zooplankton, the number of sequences in each OTU was taken as the abundance of OTUs (Zhan et al., 2014).

| Statistical analyses
All the following analyses were performed using R (R Core Team, 2018).
Three indices of α-diversity were calculated, namely Shannon index, inversed Simpson index and observed number of OTUs. The variation of community structure in five habitats was compared using analysis of similarity (ANOSIM) and illustrated by non-metric multidimensional scaling (NMDS) ordination based on Bray-Curtis distance. The variation of all environmental variables in five habitats was also compared using ANOSIM and illustrated by NMDS based on Euclidean distance.
Prior to statistical analyses, environmental variables, except pH, were log 10 (x + 1) transformed to improve normality and homoscedasticity. To explore the correlation between community dynamics and environmental variables, Mantel test was conducted to select variables which were significantly correlated with community variation. Redundancy analysis (RDA) was conducted to further explore the influence of selected environmental variables on community composition. In addition, forward selection was performed to select relatively more important variables to build a parsimonious RDA model (Blanchet, Legendre, & Borcard, 2008). Spearman's rank correlation analysis was conducted to further verify the influence of the crucial variables on OTUs which exhibited a high correlation with the first two axes of RDA (i.e., high-fitness OTUs).

| Construction and analysis of cooccurrence networks
To illuminate the relationships between bacterioplankton, zooplankton and environmental variables, a full phylogenetic molecular ecological network (pMEN) was constructed (Deng et al., 2012;Zhou et al., 2010;Zhou, Deng, Luo, He, & Yang, 2011). This method uses the random matrix theory (RMT) to identify the appropriate similarity threshold automatically before network construction. For both communities, only OTUs appearing in >50% of the total samples were used for network computation (Deng et al., 2012;Zhou et al., 2010). The relative abundance of OTUs and environmental variables was log 10transformed, and missing values were filled with a very small number (0.01) if paired valid values were available. A symmetric correlation matrix was calculated by the Spearman's correlation coefficient. Only the similarity values above a certain threshold were remained for calculating matrix eigenvalues. Random network was generated by a randomization procedure with 100 permutations to evaluate whether or not the constructed network was random (Deng et al., 2012). In the network, positive associations indicate interspecies cross-feeding, coaggregation or niche overlapping, whereas negative associations suggest competition, or niche portioning and/or resistance to be grazed (Bartomeus et al., 2016;Faust & Raes, 2012). The construction and analysis of network were performed using the pipeline MENA (http:// ieg4.rccc.ou.edu/mena; Deng et al., 2012). The statistical results of associations (edges) between bacterioplankton and zooplankton were visualized with the online version of Circos (Krzywinski et al., 2009).

| Alpha diversity and taxonomic composition of plankton communities
High-throughput sequencing of bacterioplankton community yielded 1,362,628 high-quality sequences, resulting in 7,403 OTUs.
Rarefaction curves indicated that current sequencing depth was sufficient to capture most biodiversity ( Figure S1a). No significant variation was observed for three α-diversity indices across all five habitats (Figure 2a-c Geographical distribution of the observed phyla or classes was uneven in five habitats (Figure 3a). For example, the abundance of Betaproteobacteria was higher in PW and RW2 than that in other habitats, while the abundance of Bacteroidetes was higher in MW, LW and PW than that in other habitats. The abundance of Planctomycetes was highest in MW, and the abundance of Cyanobacteria was lowest in PW. Geographical distribution of lower rank taxa was also uneven in five locations, such as at the family level ( Figure S2a). For example, the abundance of Burkholderiaceae (Betaproteobacteria) was higher in PW and RW2 when compared to the other habitats, while the abundance of Chitinophagaceae (Bacteroidetes) was higher in MW and PW than that in the other habitats.
High-throughput sequencing yielded 1,109,762 high-quality sequences, resulting in 438 OTUs for metazoan zooplankton. Most samples reached the saturation stage, implying that current sequencing depth was sufficient to capture most biodiversity ( Figure S1b).
Similarly to bacterioplankton, no significant variation was observed for three α-diversity indices across all five habitats (Figure 2d-f). At F I G U R E 2 Alpha diversity of bacterioplankton (a-c) and zooplankton (d-f) communities. LW, lacustrine wetland; MW, marsh wetland; PW, paddy wetland; RW, riparian wetland all sites, the dominant phyla of zooplankton community were Rotifera (45.52%) and Arthropoda (44.22%). The distribution of zooplankton phyla was uneven in five habitats (Figure 3b). For example, the abundance of Rotifera was higher in MW and RW1 than that in the other habitats, while the abundance of Arthropoda was higher in LW and PW than that in the other habitats. Distribution at lower taxa rank was also uneven in five habitats, such as at the class level ( Figure S2b).

| Influence of environmental variables on community diversity
Similarly to plankton communities, environmental variables in five habitats were significantly different ( Figure S4 and Table S2). For example, the concentration of COD, TN and other nutrients was lowest in LW, implying relatively better water quality in this habitat (Table S2). Mantel test indicated that 15 out of 25 environmental variables were significantly correlated with bacterioplankton community. Forward selection further revealed pH was the most important variable that structured community composition of bacterioplankton, followed by As, TDS, Fe, Mn and Na (Table 1). A parsimonious RDA model was built with these six variables (p = 0.001; Figure 5a). These selected variables explained 31.4% of the total variation at the OTU level. Correlation analysis indicated that these variables exerted distinct influence on different taxa (Table S3).
Mantel test revealed 14 out of 25 environmental variables were significantly correlated with the composition of zooplankton communities. Forward selection further revealed NH 4 as the most important variable to structure zooplankton communities, followed by As, TDS, Fe, Mn and Na (Table 1). A parsimonious RDA model was built on these variables (p = 0.001; Figure 5b). These five variables explained 29.8% of total community variation. These variables exerted discrepant influence on different taxa (Table S4).

| Complex interactions depicted by cooccurrence networks
To investigate the association among bacterioplankton, zooplankton and environmental variables, a scale-free network (R 2 of power law: 0.952) was constructed. Multiple parameters showed that the empirical network was significantly different from relevant random network, suggesting that the observed interactions were non-random (Table 2). F I G U R E 3 Phylum level composition of bacterioplankton (a) and zooplankton communities (b). LW, lacustrine wetland; MW, marsh wetland; PW, paddy wetland; RW, riparian wetland In total, the network was consisted of 373 nodes, including 301 bacterioplankton units, 60 zooplankton units and 12 environmental units.
We identified 646 taxon-taxon edges and only 32 environment-associated edges, suggesting that a majority of community variation should be explained by biotic interactions. Among environmental factors, the frequent drivers of network connections were Na, As, pH and TP. Among the taxon-taxon interactions, positive associations (85%) outnumbered the negative exclusions (15%). Besides, we observed a non-random edge distribution with regard to phylogeny ( Figure 6).
For example, most positive interactions were observed within certain groups, especially for Rotifera and Arthropoda in zooplankton. Positive associations were enriched in Actinobacteria, Alphaproteobacteria and Planctomycetes, which were dominant components of bacterioplankton communities. Although negative associations were also enriched in the three taxa, they were observed mostly between different taxonomic groups. The cross-kingdom association (21/646) between bacterioplankton and zooplankton was less than within-kingdom association (593/646).

| D ISCUSS I ON
In order to facilitate the conservation and restoration of wetland biodiversity, it is crucial to deeply investigate potential mechanisms that structure biodiversity in wetlands (Erwin, 2009;Tockner & Stanford, 2002;Xu et al., 2019;Zedler, 2000). In this study, we observed distinct community composition and structure of both bacterioplankton and zooplankton in habitats with varied levels of anthropogenic disturbance (i.e., protection status). The environmental filtering process played an important role in driving the dissimilarity of communities, with environmental variables explaining approximately 30% of community variation (Table 1). Interestingly, the two types of communities were influenced by different factors, suggesting environmental stressors in the same habitat could distinctly influenced different communities ( Figure 5; Table 1). The interactome-based analyses added extra explanation for the community variation. Therefore, this study successfully dissected ecological processes and mechanisms which drove the community structure and geographical distributions of plankton biodiversity, and the findings in this study are expected to benefit biodiversity conservation and restoration in wetlands.

| Variability of plankton communities in different habitats
Characterizing the composition and geographical distribution of communities is fundamental to assess degradation status and to establish rationale targets for wetland biodiversity restoration (Fang et al., 2006;Meli et al., 2014;Zedler & Kercher, 2005). Wetlands naturally present a high level of habitat heterogeneity, and as a result, most local populations of wetland species are small and isolated and thus vulnerable to extinction (Gibbs, 2000). Numerous pioneering studies have revealed that species, especially those flagship species such as waterbirds, were distributed in a few suitable habitats after anthropogenic disturbance (i.e., habitat heterogeneity derived from human activities; Fang et al., 2006;Yang et al., 2016;Wang et al., 2017). However, relatively less attention was paid to the composition and geographical distribution of microscopic communities. Due to the relatively large population size and high dispersal potentials, microscopic organisms have long been considered to be randomly   (Lindström & Langenheder, 2012;Martiny et al., 2006).
However, a large number of recent studies have repudiated this old viewpoint by revealing the non-random geographical distribution of these communities (Lindström & Langenheder, 2012;Soininen et al., 2007;Thompson et al., 2017), even in running water ecosystems where significant environmental gradients exist (Peng, Xiong, & Zhan, 2018;Xiong et al., 2017;Yang et al., 2018). Similarly, we also detected significant geographical variability of bacterioplankton and zooplankton communities in Sanjiang wetland in this study (p < 0.05; Figure 4). The results indicate that, similarly to many macro-organisms, the microscopic organisms can be largely threatened by anthropogenic disturbances, as environmental factor-induced species sorting promotes the geographical isolation and prohibits possible rescues from neighbouring populations.
For both types of communities in Sanjiang wetland, the species richness (i.e., α-diversity) did not significantly vary among habitats ( Figure 2), while the species composition and relative abundance (i.e., β-diversity) varied greatly (Figure 3). The provision of ecosystem services is supported not just by species richness, but also by species abundance and geographical distribution (Meli et al., 2014). Therefore, a comprehensive assessment of different biodiversity components is recommended to reflect the succession dynamics of wetland biodiversity. Regarding community composition at the OTU level, the percentage of unique OTUs in each habitat was smallest in lacustrine wetland (LW) habitat for both bacterioplankton and zooplankton communities ( Figure S3). The lower level of uniqueness of OTUs in LW suggests a higher level of migration to surrounding areas and thus LW has the potential to serve as crucial seed banks for biodiversity restoration. This result also indicates that nature reserves are crucial to maintain biodiversity, since LW habitat sits at the core of Sanjiang Wetland National Nature Reserve and is under better protection than the F I G U R E 5 The ordination plot of redundancy analysis (RDA) of bacterioplankton (a) and zooplankton (b) communities. The coloured circles represent sites that were enclosed by polygons. Operational taxonomic units (OTUs) with high fitness (>70% for bacterioplankton and >50% for zooplankton) were displayed to improve clarity. Red arrows indicated environmental variables, while black arrows indicated OTUs. LW, lacustrine wetland; MW, marsh wetland; PW, paddy wetland; RW, riparian wetland other habitats. Significant variation was also observed at the taxonomic level, such as the reverse relationship between Rotifera and Arthropoda (mostly crustaceans; Figure 3b). Such a reverse relationship was common in aquatic ecosystems (Xiong et al., 2017;Yang et al., 2018) and was mainly related to interspecific interactions, such as competition or predation (Laxson, McPhedran, Makarewicz, Telesh, & MacIsaac, 2003;Meyer, Hampton, Ozersky, Rusanovskaya, & Woo, 2017).

| The importance of environmental filtering
Studies have summarized that abiotic factors impact community composition and geographical distribution mainly through two ecological processes-environmental filtering and dispersal limitation (Heino, Melo, & Bini, 2015;Leibold et al., 2004;Soininen et al., 2007). In this study, several environmental variables were identified to drive the distinct community composition of plankton communities in Sanjiang wetland ( Figure 5; Table 1), stressing the importance of environmental filtering. The strength of environmental filtering depends largely on the length of environmental gradient (Heino et al., 2015;Xiong et al., 2017). In Sanjiang wetland, the environmental gradient mainly comes from the significantly different water quality in five habitats, which was mainly derived from varied levels of anthropogenic disturbances. Due to long-term reclamation, large area of natural wetlands has been transformed to croplands, resulting in severe degradation and loss of wetlands (Song et al., 2014;Wang et al., 2015). Newly launched protection programmes, such as the establishment of nature reserves, have been beneficial to conservation of wetlands. For example, the relatively better water quality of LW habitat results from the better protection status (Table S2). However, the other habitats are suffering from different levels of anthropogenic disturbance, such as ongoing reclamation in riparian wetlands (RW1 and RW2 habitats), and low efficiency of restoration from croplands (PW habitat;Song et al., 2014;Wang et al., 2015). Therefore, it still remains a great challenge to restore and conserve the biodiversity of Sanjiang wetland.

| The complex interactions of biotic and abiotic factors on biodiversity
Our results showed that the two types of communities responded differently to environmental stressors in the same habitats ( Figure 5; Moreover, the taxon-taxon associations across functional types and taxonomic groups were non-randomly distributed on the network. Most of the interactions were positive (85%), suggesting that species mutualism or parasitism might be the dominant processes to regulate ecosystem functioning in Sanjiang wetland. Interestingly, positive interactions mostly occurred within each kingdom. Species with close phylogenetic relationship are supposed to share similar biological functions or occupy similar niches, and thus facilitate the uptake and metabolism of resources (Deng et al., 2012;Faust & Raes, 2012). However, negative interactions occurred mainly between bacterioplankton and zooplankton or within zooplankton. As some taxa of zooplankton could feed on bacteria or small-bodied zooplankton, grazing or predation might play important roles (Baltar et al., 2016;Fermani et al., 2013). In return, bacteria could release some secondary metabolites, such as Shiga toxin, to defend the predation of zooplankton (Jousset, 2012;Matz et al., 2008). Therefore, the complex interactions, especially taxon-taxon associations, contributed to the variation of both bacterioplankton and zooplankton communities.

| CON CLUS IONS
Effective conservation of wetland biodiversity is hampered by the lack of a comprehensive understanding of ecological processes and mechanisms that influence biodiversity. In this study, we found complex interactive effects of abiotic and biotic factors to drive the distinct distribution patterns of plankton communities in Sanjiang wetland. The selective strength of abiotic variables was mainly derived from different levels of anthropogenic disturbances.
Interestingly, bacterioplankton and zooplankton communities responded differently to variables in the same habitats, suggesting that a comprehensive survey is needed to take diverse organisms into consideration when restoration programmes are established.
Furthermore, the incompleteness of abiotic variables to predict community dynamics was partially supplemented by the integration of biotic interactions, calling for attention towards the interactomebased framework to examine the synthesizing impacts of both abiotic and biotic variables. Outcomes from comprehensive surveys can help identify conservation priorities and improve the predictability of conservation programmes.

DATA ACCE SS I B I LIT Y
These sequence data have been submitted to the GenBank Sequence Read Archive under the accession numbers PRJNA 542032 (zooplankton) and PRJNA 542034 (bacterioplankton).

Aibin Zhan
https://orcid.org/0000-0003-1416-1238 F I G U R E 6 Taxonomic patterns within the co-occurrence network. Top nine interacting taxon groups were depicted as coloured segments in a Circos plot, in which ribbons connecting two segments indicate (a) co-presence and (b) exclusion links. Size of the ribbon was proportional to the number of links (co-presence and exclusions) between the operational taxonomic units (OTUs) assigned to the respective segments. Links were dominated by Actinobacteria, Alphaproteobacteria and Planctomycetes