Microbial co-occurrence networks as a biomonitoring tool for aquatic environments: a review

Healthy aquatic ecosystems are pivotal for supporting life on Earth, and of immense social, economic, cultural and environmental value. However, the health of aquatic ecosystems is being put at risk, due to the direct and indirect effects of human activities (Halpern et al. 2008). With the growth of the human population and the expansion of cities globally, pollution, alterations to land use, climate change, overexploitation and invasive species are having an increasingly pronounced detrimental effect on aquatic systems (Friedman et al. 2020). Individually, and in concert, these anthropogenic stressors are placing unprecedented pressure on all aquatic ecosystems, eroding biodiversity and altering ecosystem functions and services at an unparalleled scale and pace (Steffen et al. 2011). Consequently, a deeper understanding of the potential impacts and extent of human activities on aquatic ecosystems is pertinently needed to better inform management actions that are capable of both protecting essential water resources, while also supporting the growth and expansion of the human population.

Freshwater, estuarine and marine ecosystems are complex, heterogeneous, and dynamic environments (Heino et al. 2015). They are host to a great variety of species, with the majority of their biomass, up to 90% (Costello et al. 2010; Snelgrove 2010) being composed of microorganisms. Microorganisms, most notably prokaryotes, provide invaluable functions including carbon and nitrogen fixation and the remineralisation of organic matter, as well as forming the basis of oceanic food webs (Jeffries et al. 2016). The global carbon and nutrient cycles are thus dependent on the activity of these microorganisms (Azam and Malfatti 2007). As such, human induced changes to aquatic microbiota can have secondary effects on productivity, food-web structure and energy flow, and carbon export (Sagova-Mareckova et al. 2021). Moreover, the microbiota of these systems are highly responsive to local and global pressures, because of their high metabolic and growth rates (Allison and Martiny 2008). Consequently, the microbiota can be highly influenced by the input of both organic and inorganic pollutants (Holman et al. 2021; Sagova-Mareckova et al. 2021). For example, it has been found that in an aquatic mesocosm, bacterial community composition, structure and function can be markedly altered by increasing concentrations of copper (Sutcliffe et al. 2019; Codello 2021).

Microbial responses to stress, such as changes diversity, composition and structure can be used to detect disturbances in water ecosystems (Williams et al. 2014). Although the use of aquatic biota as indicators of disturbance and ecological condition is an integral part of water quality and sediment monitoring (Sun et al. 2012), to date, they are very eukaryotic centric, with the routine monitoring of microorganisms, including bacteria, archaea and protists, rarely employed (Bohan et al. 2017; Sagova-Mareckova et al. 2021).

Indicators based on microorganisms can be selected from a broad assortment of taxonomic groups and cover a wide variety of functions and services (Aylagas et al. 2016). Microorganisms are highly sensitive and have an overall population that include both niche specialists and generalist microbes (Bell and Bell 2021). The census of aquatic microorganisms is today possible due to the development of molecular technologies like high throughput sequencing and readily available genetic databases. These advances, which have allowed the construction of effective new approaches such as the culture-free analysis of a microorganism’s genetic material (Shendure and Ji 2008) enable numerous samples of complex assemblages to be sequenced simultaneously.

Analysis of microbial genetic material with metabarcoding and amplicon sequencing techniques is accurate, fast, cost-effective and capable of analysing thousands of environmental samples at the same time, providing a better understanding of microbial communities than previous techniques. These methods can reveal information on the most abundant microbial species in aquatic environments, on rare species, recently introduced invasive species, and can report the presence of organisms that cannot otherwise be sampled or identified with more traditional techniques (Aylagas et al. 2016; Cristescu and Hebert 2018).

To date, the focus has been on composition and univariate endpoints, which does not inform how species interact with each other and as a collective. Most microbial communities are complex systems with intricate interaction patterns showing a collective behavior that cannot be identified by just looking at the single entities or composition (Cimini et al. 2019). Consequently, there is a need to use more advanced analytical methods such as co-occurrence network analysis and their properties to monitor the biological environment and provide a far more complete picture (Tylianakis et al. 2010).

In this review, we demonstrate the potential utility of the microbial community co-occurrence networks and their properties as a biomonitoring tool to assess aquatic health. To date, a number of studies has employed this approach to characterise the system-level responses to environmental change in aquatic systems, including anthropogenic disturbances (Bohan et al. 2017; Cordier et al. 2020; Di Battista et al. 2020). However, it has yet to be used to its full potential because of a lack of empirical data and conflicting views with regards to the ecological implications of any perceived changes in network structure and its associated metrics (Karimi et al. 2017; Derocles et al. 2018; Barroso-Bergada et al. 2020). To address this, firstly, we provide a summary of aquatic microbial networks, specifically co-occurrence networks, and their main structural components, such as nodes and edges. Secondly, we describe network properties, and where possible, what ecological information can be drawn from these networks. Finally, we illustrate the use of co-occurrence network as a biomonitoring tool by providing an overview of existing studies, while articulating the potential advantages and pitfalls associated with applying co-occurrence networks for biomonitoring in aquatic systems.

The theory that underpins ecological co-occurrence analysis originates from Jared Diamond’s seminal work on the assembly rules of ecological communities, which predicted that competitive interactions between bird communities in New Guinea would result in non-random co-occurrence patterns (Diamond 1975). The bird population was found to follow a ‘checkerboard’ distribution representing one of the first pieces of evidence that competition can play a key role in determining which species are found within a community. A network is a representation of a community consisting of a multitude of parts, also called nodes, which are connected by a binary interaction that allows exchanges or communications between the parts (Newman 2010). This concept has been used by many fields with myriad applications, ranging from social networks (e.g. Facebook), the World Wide Web to electrical grids. Some of the most common types of biological networks are protein–protein interaction networks, metabolic networks, genetic interaction networks, gene and transcriptional regulatory networks, cell signaling networks, food webs, and networks of interactions between species (co-occurrence) (Newman 2010). Although the biological applications are broad and capture a wide range of systems and processes, they often have a common characteristic, being scale-free networks, also called real-world networks, which inherently contain a relatively low number of nodes that are highly connected to other nodes (Fig. 1) (Bianconi and Barabási 2011).

Fig. 1.  Graphical representation of a scale-free network. Nodes are represented by green dots and edges are represented by grey lines between the nodes. The graph was plotted with Mathematica (ver. 12.1, see https://www.wolfram.com/mathematica/) using the function Barabasi Albert Graph Distribution (Wolfram 1999).

Biological networks are generally compared to a random network, where it is rare to find nodes that have significantly more or fewer connections than the average node (Barabasi and Bonabeau 2003). Most biological complex systems have a similar architecture and similar organisational principles. Consequently, the scale-free structure has important implications, for example, in controlling the spread of viral diseases such as HIV and COVID-19 (Barabasi and Bonabeau 2003; Murali et al. 2011; Pal et al. 2020). Similarly, the interaction between microorganisms have been shown to have scale-free distribution in many freshwater and marine environments (Ruan et al. 2006; Steele et al. 2011; Kara et al. 2013).

Although there are numerous types of networks, here we focus on co-occurrence networks. Co-occurrence networks are derived from the presence or abundance and interactions of taxa occurrences within a repeatedly sampled unit e.g. multiple samples from a site. Co-occurrence networks are founded on the interactions between taxa derived from significant bipartite relationships. The two fundamental components are nodes and edges. Nodes represent taxa that have a significant interaction with other taxa, and edges represent the co-occurrences bipartite interactions between taxa, allowing for the possibility to represent and study the interactions of the whole microbial community (Fig. 2). A microbial network co-occurrence analysis allows a holistic study of the environment that not only focuses on presence or absence and composition of organisms in the ecosystem but also on how they interact (Layeghifard et al. 2018).

Fig. 2.  A graphical representation of an aquatic sediment microbiome co- occurrence network (Codello 2021). Nodes are represented by green dots, and edges with grey lines. The network was calculated with SpiecEasi and represented with iGraph.

A holistic methodology can facilitate our understanding of the complex interactions that occur within biological systems. This is evident in the case of marine bacterioplankton, where ten different ‘ecological species’ of SAR11, phylogenetically similar, have been found due to the unique co-occurrence network relationship between them, other organisms and environmental parameters (Fuhrman and Steele 2008). For a global biomonitoring approach, networks can provide a comprehensive understanding of how ecosystems function across systems and biogeographical regions (Raffaelli 2007; Bohan et al. 2017). Thus, it has been suggested that inclusion of information that examines the interactions between species can greatly benefit our understanding of community ecology, and how communities respond to change (Valiente-Banuet et al. 2015).

An interaction between bacterial species involves three main exchanges: food transfer, information transfer, and gene transfer (Fuhrman 2009). Food transfer in microorganisms creates a food web that is at the centre of the marine ecosystem, allowing the exchange of energy and biochemical compounds (Steele et al. 2011). Information transfer also occurs in microbes since they can receive information cues from the external environment (such as a chemical signal from other microbes) and are able to respond (Westerhoff et al. 2014). Gene transfer can happen both in the form of gamete exchange and lateral gene transfer, e.g. resistance cassettes and antibiotic resistance genes (Gillings et al. 2008; Fuhrman 2009).

Nodes represent the first fundamental component of a network. Nodes, are the objects at the start and ends of the connection, and are commonly represented as Operational Taxonomic Units (OTUs), Amplicon Sequence Variants (ASVs), taxa, species, genera, families, classes, phyla or kingdoms, or functional groups (Karimi et al. 2017). The role of the node in relation to the rest of the network allows us to understand the dynamics of the system. ‘Hubs’, ‘bottle necks’ and ‘keystone’ species represent the three main node roles. A ‘hub’ defines a node that has a high level of degree centrality (a high number of connecting neighbouring nodes) (Paine 1995; Delmas et al. 2019), whereas a ‘bottleneck’ is a node that has a fundamental role in connecting other nodes (Peura et al. 2015). Usually, ‘bottlenecks’ are between nodes and get crossed by many pathways inside the network; they have a high betweenness centrality and high transitivity (Peura et al. 2015; Delmas et al. 2019). Common marine and freshwater taxa are typically ‘hub’ and or ‘bottleneck’ species likely due to their ability to adapt to environmental change (Peura et al. 2015; Lin et al. 2019). Species that have a large effect on the community are called ‘keystone’ species, they are less abundant than other organisms in the community, but the other nodes depend on their effect (Power and Scott Mills 1995; Berry and Widder 2014). These keystone species play a critical role in the stability of the system (Peura et al. 2015). Apex predators are an example of a keystone species in macro-ecology, their predatory behaviour assists in the control of the population size of their prey (Berry and Widder 2014). ‘Keystone’ species can be ‘hubs’ but ‘hubs’ are not always ‘keystone’ species. A ‘hub’ species is one that when eliminated can be substituted by another competitor species occupying the same niche with the same or similar function (Berry and Widder 2014; Delmas et al. 2019). However, when a ‘keystone’ species is eliminated, the system undergoes significant redistribution and potential collapse (Valiente-Banuet et al. 2015). It has been hypothesised that in a fragile network, the removal of 25% or less of the highly connected nodes should result in collapse, whereas in a robust and stable network the removal of half of the nodes would not change the network topology (Estrada 2007).

Collectively, rare organisms can influence important biogeochemical processes, thus, in some instances they can also be considered ‘keystone’ species. Although rare taxa may also be important, it should be emphasised that in most cases these taxa do not become nodes as they are only found in a small proportion of the samples used to create the network, and are therefore removed during the initial filtering period, or fail to have robust statistical relationships with node, they must be other taxa (Fuhrman 2009).

Environmental factors can also be incorporated in the network structure as a type of node where the edge exemplifies the effect of the factor on the microbial species node (Li et al. 2018). In seawater samples, a strong archaeal–bacterial association was found due to environmental pressure (Parada and Fuhrman 2017), whereas in thawed ponds and lakes, the importance of environmental factors was highlighted by dissolved organic carbon and conductivity being the most connected nodes (Comte et al. 2016). The inclusion of environmental factors is a powerful characteristic of network analysis that can help to describe ecological niches and discern what factors have the greatest influence over the microbiome (Chaffron et al. 2010; Gao et al. 2019; Mikhailov et al. 2019).

The second fundamental component of networks are edges, which represent interactions between nodes. These can be broadly categorised in accordance with their co-occurrence interactions (see Fig. 3): positive interactions, combination of positive or negative with neutral interactions and negative interactions (Lidicker 1979). Theoretically, positive interactions suggest commensalism and mutualism (Bronstein 2009). Commensalism describes a cooperation between microbes where only one of the two gains an advantage and the second is not affected at all (Lidicker 1979). However, mutualism describes a cooperation where both organisms benefit from the interaction. For example, the bacteria Bacteroides thetaiotaomicron is found in the human distal intestinal microbiome, where it assists in the breakdown of indigestible polysaccharides, making them bioavailable to humans (Bäckhed et al. 2005). In this example, The B. thetaiotaomicron benefits by obtaining its energy from the polysaccharides, whereas the host (human), mutually benefits as the polysaccharides are now converted to digestible form. Mutualism can also be symbiotic, i.e. when both microorganisms are required to survive (Bronstein 2009). Organisms can also have mutual relationships that are syntrophic, that is, when the metabolism of the two organisms is complementary (Barton and Northup 2011). For example, in marine gas–hydrate-rich sediments with high rates of methane-based sulfate reduction, Archaea grow in highly populated conglomerates and are surrounded by sulfate-reducing bacteria, thereby creating a syntrophic consortium that mediates methane oxidation (Boetius et al. 2000).

Fig. 3.  Schematic classification of the biotic co-occurrence interactions based on the relationship between Species 1 and Species 2. These interactions were organised along the co-occurrence patterns gradient: at the top there are positive co-occurrence interaction patterns and, at the bottom, the negative co-occurrence interaction patterns. The colour of the dot represents if a species benefitted (green), was harmed (red) or was not affected (blue) by the interaction with the other species. Green arrows represent positive interactions, red arrows represent negative interactions.

Positive–negative interactions occur when one taxon gains an advantage and the partner is negatively affected, e.g. predation and parasitism (Barton and Northup 2011). In predator–prey interactions, the predator improves its condition by altering the condition of the prey (Lidicker 1979). For example, in marine sediments, Pseudoalteromonas was found to prey on Gram-positive bacteria by secreting pseudoalterin. Pseudoalterin binds to and degrades the glycan strands of Gram-positive bacterial peptidoglycan. This causes the cell to die and release nutrients used for the growth of the predator (Tang et al. 2020). Parasitism represents a host–parasite interaction where the parasite is physically associated with the host and it is in the privileged position to gain an advantage at the expense of the host (Lidicker 1979). In marine environments, abundant viruses can infect cyanobacteria such as Synechococcus sp. Electron microscopy analysis showed that 1–3% of the native cyanobacteria contained viruses close to the end of the lytic cycle. It seems that cyanophages are probably responsible for the death of up to 3% day⁻¹ of Synechococcus in temperate Atlantic waters (Fuhrman 1999). Finally, negative interactions suggest amensalism or competition (Barton and Northup 2011).

Amensalism is characterised by an interaction where one organism is affected by the presence of the other, but the second does not gain any advantage (García et al. 2017). This type of interaction has been observed in the fermentation of residual cheese whey by the Pseudomonas taetrolens and Lactobacillus casei collaboration. During their interaction there was no variation in Lactobacillus, whereas the Pseudomonas community had reduced growth (García et al. 2017). Within the co-occurrence context, competition is an interaction where both nodes are adversely affected, for example, when they directly compete for the same resources (Barton and Northup 2011). An interesting example of bacterial competition was studied between Pseudomonas sp. and Scenedesmus quadricauda, which competed for inorganic phosphorus under laboratory conditions. The dynamics between the two bacteria suggested that, at a steady state with sufficient carbon, Pseudomonas sp. reduced phosphorus below the level required by S. quadricauda, whereas when Pseudomonas sp. did not have enough available carbon S. quadricauda prevailed (Grover 2000).

Although the interpretations of statistically derived interactions observed in co-occurrence networks are driven by ecological theory, in reality, biological relationships between interactions and co-occurrences are often ambiguous (Blanchet et al. 2020). This is particularly true for aquatic microbial communities, where some species are naturally rare or transient, and therefore only interact under specific situations (Blanchet et al. 2020). Furthermore, macro-ecological networks have shown that co-occurrences reflect more niche preferences rather than biotic interactions (Freilich et al. 2018). It is safe to assume that bacteria can reside in the same habitat without interacting so, inference based on network associations, in some instances, may result in erroneous results (Röttjers and Faust 2018). Another issue is that, although the gradient of diversity is continuous, the categories we are describing are discrete. Specifically, it is challenging to discern between parasitism and predation, as well as understanding the difference between mutualism and commensalism, or between amensalism and competition. This is because it is difficult to capture the positive gradient between mutualism and commensalism and the negative gradient between amensalism and competition (Röttjers and Faust 2018). Originally, most networks were inferred from presence–absence data; however, now abundance (count) data provide more information to refine the inferences (Blanchet et al. 2020).

In networks based on metabarcoding data, edges represent the correlation value between OTUs or ASVs. The network edges can be calculated using a wide range of methods from simple correlations that are quick to calculate, to more complex approaches, which require higher computational power. The simplest methods are based on dissimilarities, such as Bray–Curtis and Kullback–Leibler dissimilarity indices, and correlations, such as Pearson and Spearman correlation coefficients (Layeghifard et al. 2017). Multiple regression and probabilistic-based methods, such as Bayesian networks or Markov networks, can also be used to capture more complex interactions (Layeghifard et al. 2017). One of the great challenges, especially when using NGS-derived bacterial data, is addressing biases associated with composition. There are a number of packages specifically designed to minimise the artefacts of such data, including SparCC, MIC, SPIEC-EASI, CoNet and others (Faust and Raes 2012; Layeghifard et al. 2017). An additional important component of network analysis is the capacity to visualise the data. Again, there are numerous approaches, including Cytoscape (Shannon et al. 2003) and iGraph (Csardi and Nepusz 2006), as well other functions within R. The choice of approaches used to produce and analyse a network may influence the overall metrics produced by a network and subsequent interpretation. This is highlighted in Codello (2021), where the same system was examined using SPIEC-EASI and CoNet. Although the findings were not conflicting, they highlighted that the metrics produced from the two approaches varied in their sensitivity to detect stress.

In the aquatic environment, environmental factors play a key role in shaping microbial communities (Mikhailov et al. 2019). These include water stratification (Zaccone and Caruso 2019), different levels of light penetration (Parada and Fuhrman 2017), temperature gradient, seasonal turnover (Gilbert et al. 2012; Mikhailov et al. 2019), spatial heterogeneity (Fuhrman 2009), and disturbances present in the water (Liu et al. 2015), which can change the structure or topology of the network. For example, in eutrophic freshwater Lake Taihu, water temperature, chlorophyll a and phytoplankton density were the main factors driving the variation of bacterial community composition and structure between seasons (Zhu et al. 2019). Co-occurrence networks formed larger networks with more nodes in spring and summer than in autumn and winter. Other studies confirmed that warmer seasons support the activity of more bacteria than colder months in water systems (Crump and Hobbie 2005; Gilbert et al. 2012; Kara et al. 2013). However, after the summer network, the autumn network had the largest number of edges, suggesting that the microbial community transition from warmer temperature to colder temperature stimulates the creation of new interactions, be it positive or negative (Zhu et al. 2019). Consequently, the overall structure and topology of networks may provide insight into how a community changes.

Environmental networks are modular but contain a hierarchical structure. At the lowest level are individual nodes. A group of nodes can form small subgraphs called motifs, together motifs create modules and finally modules form a complete network (Mason and Verwoerd 2008). It has been suggested that this hierarchical structure could be the result of co-evolution, natural selection, habitat heterogeneity, the specificity of interaction and phylogenetic relatedness (Deng et al. 2012). For example, Liao et al. (2019) found co-occurrence networks derived from sediments within the tidal zone produced distinctive modular structures that could visually differentiate them from non-tidal and inter-tidal networks. Hence it has been suggested that the structural modularity of a network may also reflect the stability and resilience of the system (Pavlopoulos et al. 2011; Deng et al. 2012; Delmas et al. 2019; Di Battista et al. 2020). In planktonic prokaryotic community profiles from the Taihu Watershed in China, the co-occurrence network was grouped into four major ecological modules where the OTU nodes had a robust co-occurrence correlation (Liu et al. 2020a). The different modules performed various functions, for example, the first module was attributed to have four functions: methanotrophy, aerobic ammonia oxidation, nitrification and ureolysis; the second module had functional groups related to metabolism and human health; and the third module had an association with chlorate reduction, methanol oxidation and methylotrophy. Hereby, clearly demonstrating that in this case, functional groups were linked by ecological associations (Liu et al. 2020a). The analysis of microbial networks can be coupled with specific functions, which allows us to have a general understanding of the microbial phylogenetic distribution, their correlations, and their collective functions (Bohan et al. 2017). The analysis can be integrated with genetic, biochemical, and chemical functions (Bowen et al. 2013). Genes with known functions can profoundly influence the structure of networks, and consequently, in some incidents where the function of a gene is unknown, its role might be revealed by the taxa it interacts with (Galand et al. 2018). For example, if a taxon is highly connected with a group of nitrogen-fixing microorganisms, there is a high probability that it too will be a nitrogen fixer, thus creating a direct link between the taxon’s role in a network and their functional role (Bohan et al. 2017). Hence, determining the network structure of co-occurrence interactions can help define ecological niches and characterise microbial ecotypes.

Mathematical approaches are not only used for network graphical representations, but to understand the behaviour of networks, enabling the properties of a network to be summarised using a range of metrics (Proulx et al. 2005). Network properties are features that can provide information on the biological system at the base of the network. There are specific properties to describe and measure network topologies, and can be used to reveal patterns, and inform the structure of the microbial co-occurrence network within aquatic ecosystems (Steele et al. 2011; Röttjers and Faust 2018). The subsequent section provides an overview of some of the key network properties, and how these properties specifically relate to microbial communities, and where possible, in response to disturbance. The basic properties that can describe the whole network and have potential to be applied within a biomonitoring framework include: the number of nodes, number of edges, number of components, diameter or radius, degree distribution, average path length, density, cluster coefficient and centrality measures.

The number of nodes and number of edges measures the total number of nodes and the total number of connections between nodes respectively. In the lake sediment bacterial community during a toxic bloom event caused by untreated sewage, algal bloom and aquaculture, the network structure increased its dimensions. Before the bloom, the network had 100 nodes and 1851 edges, compared to 200 nodes and 3193 edges after the bloom. This expansion highlighted the interference of the toxic bloom on the bacterial community composition (Zhou et al. 2021). It is possible that at first the bacterial community formed a structure that was disrupted by the toxic bloom. The toxic event and the proliferation of the algal bloom increased the microbial complexity with more nodes and increased the number of negative interactions between the bacterial community and the algal bloom (Zhou et al. 2021).

The number of components is the number of modules or groups of nodes that are connected between each other but not to other nodes of the same network, i.e. there are no possible paths between nodes from different components (Newman 2010). In essence, it is a measure of completeness. When there is the fragmentation, former relatively larger networks become dissociated into smaller disconnected networks, resulting in an increase in the number of components (Pavlopoulos et al. 2011). This phenomenon was seen in a river system where the number of components increased downstream compared to upstream because of hydrological changes and metacommunity dynamics (Widder et al. 2014). In this example, it was hypothesised that the elevated physical disturbance related to flow patterns and sedimentary dynamics in the downstream river, caused the loss of keystone taxa, thereby increasing fragmentation in downstream sites (Widder et al. 2014).

The diameter or radius is a measure of the dimension of the network and it is the longest of all the path lengths. It measures how quickly a community can respond to a stimulus (Pavlopoulos et al. 2011). The network constructed with river bacterial communities data changed diameter length according to seasons: spring had the shortest diameter, describing a very efficient and quick to respond community; and autumn had the longest diameter, indicating a slowing down of community collaborations (Lin et al. 2019). These changes suggest that seasonality interferes with the structure of the bacterial network community and with the efficiency to collaborate (Lin et al. 2019).

The average path length (N_ij) (Eqn 1) of a pair of nodes is the average number of edges (E) that have to be crossed between two nodes (i, j) to have the shortest distance (d_ij) (Newman 2010).

The average path length is another measure of the expediency of the network reaction to a change in the environment or to a change in the node composition (Newman 2010). The length informs how fast two nodes connect and pass information between them, a short path means that the system is ready to react (Mason and Verwoerd 2008). In network experiments from five aquatic time series in dystrophic and eutrophic lakes and open ocean sites, the average path lengths between nodes were shorter in empirical than in random networks (Peura et al. 2015; Liu et al. 2020b). This suggests that in aquatic environments, the average path length maybe comparatively shorter than the whole network dimension (Newman 2010).

Density (D) (Eqn 2) describes how compact or crowded the network is in relation to the number of connections per nodes, and thereby is a measure of complexity (Landi et al. 2018).

where E is the number of edges and n is the number of nodes (Newman 2010).

A sparse network has fewer interactions, lower complexity and needs less energy to be maintained, but a dense network has more interactions, higher complexity and higher energy requirements (Pavlopoulos et al. 2011). Experimental data has shown that biological networks have a lower density than random networks. In the ‘Domino hypothesis’ a less dense network is considered to be more resilient and efficient, with lower costs to maintain the system than a dense network (Leclerc 2008). By contrast, the ‘Insurance hypothesis’ suggests that a denser network has more redundant nodes, which could make the network more resistant to change (Yachi and Loreau 1999). In benthic microbial networks produced from affected offshore oil and gas drilling sites, the density was higher in non-affected sites than those closer to the drilling operations (Laroche et al. 2018). This suggests that there was a reduced complexity in the microbial co-occurrence network when the environment was disturbed by human activities and small spills of oils (Laroche et al. 2018). These two hypotheses highlight the contrary viewpoints when extrapolating the ecological relevance of variations in some network properties.

The cluster coefficient (Cl) (Eqn 3) represents the predisposition to form clusters or modules (Newman 2010). In a graph, a cluster is formed when one or more nodes have a high number of edges E and a low node degree k.

Clusters are common in natural environment networks more than in random networks (Newman 2010). There are network subsets in which species frequently interact between themselves; however, the subset interacts infrequently with species outside of the cluster (Pavlopoulos et al. 2011). A study in agricultural soils, where a single application of alkaline stabilised biosolid (used to destroy pathogens) was applied 10 years prior to sampling, found a 39% increase in the clustering coefficient, in comparison to the control soils (Price et al. 2021). This suggested that the agricultural additive (alkaline-stabilised biosolid) had a strong legacy effect on the soils and on the network structure increasing its complexity.

Centrality measures are properties that can quantitatively describe microbial interactions and identify the most important taxa and offers suggestion on the role or function of the taxa in relation to the rest of the network (Zamkovaya et al. 2021). They are: degree centrality, closeness centrality, betweenness centrality, eigenvector centrality and transitivity (Girvan and Newman 2002). They can be studied as mean values at the level of the full network, or singularly for each node.

The node degree or degree centrality (k_i) (Eqn 4) is the number of edges or connections that a node (i) has with other nodes, where n is the total number of nodes (Newman 2010)

A node with high degree has many connections to other nodes, with a lower degree, inferring a sparsely connected node (Pavlopoulos et al. 2011). Owing to the characteristic distribution of scale-free networks, most microbial communities have many nodes with few links and few nodes having a large number of links, also called ‘hubs’, making the network possibly more resistant to random failure, but also vulnerable to organised attacks to the key nodes (Barabasi and Bonabeau 2003). An example of degree centrality in aquatic sediment environment is represented in Fig. 4.

Fig. 4.  Representation of degree centrality in co-occurrence network from aquatic sediment microbiome calculated with SpiecEasi and plotted with iGraph (Codello 2021). The red colour gradient represent the degree centrality measure (from low degree centrality in white to high degree centrality in dark red).

Closeness centrality (C_i) (Eqn 5) describes how close a node is to other nodes. It is defined as:

where n is the total number of nodes and d_ij the shortest distance between i and j (Pavlopoulos et al. 2011). A node has a high closeness when it can communicate quickly with another node (Newman 2010). A key player in a network has high closeness centrality, because it has to quickly interact with many network components (Pavlopoulos et al. 2011). An example of closeness centrality in an aquatic sediment network is illustrated in Fig. 5.

Fig. 5.  Representation of closeness centrality in co-occurrence network from aquatic sediment microbiome calculated with SpiecEasi and plotted with iGraph (Codello 2021). The red colour gradient represent the closeness centrality measure (from low closeness centrality in white to high closeness centrality in dark red).

Betweenness centrality (x_i) (Eqn 6) is the degree of an individual node that lies on a path between other nodes. It is defined as:

where nⁱ_st is 1 if the node i lies on the shortest path from s to t and 0 if it does not (Pavlopoulos et al. 2011). Nodes with high betweenness connect nodes in the network that otherwise would not be connected (Newman 2010). Between centrality identifies nodes that are intermediates between neighbours and belong to many network pathways (Pavlopoulos et al. 2011). For example, species that act as a bottleneck have high betweenness. An example of betweenness centrality in aquatic sediment environment is represented in Fig. 6.

Fig. 6.  Representation of betweenness centrality in co-occurrence network from aquatic sediment microbiome calculated with SpiecEasi and plotted with iGraph (Codello 2021). The red colour gradient represent the betweenness centrality measure (from low betweenness centrality in white to high betweenness centrality in dark red).

Eigenvector centrality rank (x_i) (Eqn 7) measures the importance of the neighbouring connection. In some instances, a node’s importance is increased by the connection to other nodes that are important themselves (Newman 2010). It is defined as:

where x_i is the centrality with i nodes and the adjacent matrix A_ij (Newman 2010).

Katz centrality and Page rank are also specific types of eigenvector centrality. They are relevant for identifying key populations and a higher order of organisation (Deng et al. 2012). Not all nodes have the same importance, a high eigenvector centrality node does not have to be highly linked, and it can have less connections with more important nodes (Newman 2010).

The transitivity (T) (Eqn 8) is the probability that two nodes are connected directly or indirectly.

where l_i is the number of links of node i and m is the number of neighbours of i (Newman 2010). For example, hub nodes have high transitivity because of their many connections (Girvan and Newman 2002).

Collectively, these centrality measures are important components of network structure and integrity (Banerjee et al. 2018). For example, in the co-occurrence networks from a large eutrophic Lake Taihu (China), rare taxa, although being low in abundance, showed relatively higher values of degree centrality, closeness centrality and eigenvector centrality than abundant taxa (Zhou et al. 2021). Consequently, rare species can have important roles in the network, supporting the idea that rare species harbour genetic variability with the metabolic potential and capacity to gain a dominant role under specific conditions (Lynch and Neufeld 2015). Also, Zhang et al. (2020) identify rare bacteria as a keystone component of the community network in biochemical processes and community assemblies presenting high degree centrality and low betweenness centrality values. The keystone taxa were Denitratisoma, which played a role in the removal of nitrogen, Anaeromyxobacter, which affected the mobility of metal contaminants, and Candidatus Microthrix, which was involved in the removal of total nitrogen (Zhang et al. 2020).

Environmental change is reshaping biodiversity and the provision of ecosystem processes and services globally. The use of metabarcoding and sequencing technologies provides a new generation of biomonitoring tools, which not only capture microbial communities, but also their interactions by network-based approaches (Bohan et al. 2017; Cordier et al. 2020). Co-occurrence networks and their properties can offer new insights into ecosystem degradation and could characterise the system-level responses of environmental change, including pollution, land use, climate change, overexploitation and species invasions (Bohan et al. 2017). An added advantage of using networks is that they can be independent of taxonomy (Makiola et al. 2020), which may be particularly useful for biomonitoring in novel environments where taxonomy of biota is poorly known.

The microbial co-occurrence network is a dynamic system that changes over time and space creating or losing interactions due to abiotic and biotic pressures such as oil spills (Laroche et al. 2018), reductions in the availability of water (Hernandez et al. 2021), and longitudinal changes in rivers (Widder et al. 2014). As described above, specific network properties can indicate different responses in communities and are thus valuable for detecting and potentially revealing causality in a biomonitoring context.

The type of interactions in a network can provide information on the nature of changes in a microbial community (Liu et al. 2020c). Network properties revealed a change in soil microbiome interactions due to the water availability (Hernandez et al. 2021). Specifically, positive associations between taxa dominated the high-stress environments whereas the negative associations dominated in the lower stressed environments (Hernandez et al. 2021).

One of the most obvious changes in network structure is fragmentation. Fragmentation occurs when a single large, well-connected network degenerates to several smaller, weakly connected networks (Delmas et al. 2019). This phenomenon can be caused by anthropogenic contamination and ecological conditions. In rivers, a gradient of contamination had a gradual effect on the microbiome reflected in the disappearance of organisms (Li et al. 2018). These changes were evident in the network structure as loss of connections between microorganisms and reduction of path length. When the average path length of a network is reduced, it has been proposed that the perturbation can more rapidly spread through the network, increasing the suspectability of the network to fragmentation (Laroche et al. 2018).

Fragmentation can also occur in the network constructed from marine biofilm (Lawes et al. 2016). The biofilm has a biological structure where microorganisms stick to each other creating an extracellular matrix, building a multispecies structure to keep beneficial partners close and share functional ability (Little et al. 2008). General microbial survival rules can change when the microbes are involved in a biofilm community and they could also influence the bacterial interaction structure (Marti et al. 2013). In a study on a young marine biofilm co-occurrence network representation, composition and connectivity changed with the availability of nutrients (Lawes et al. 2016). The biofilm with ambient nutrients showed stronger connectivity, whereas the biofilm with enriched nutrients showed a fractured co-occurrence pattern (Lawes et al. 2016). These findings demonstrated a change in how the biofilm used the local resources, from a more collaborative state, where the microbiome used the energetic allocations from within the core microbiome efficiently using local resources; to a fracture network to avoid the unfavourable conditions (exposure to enriched nutrients) and did not exploit local resources to obtain energy (Lawes et al. 2016).

To be able to understand changes in the network topology it is important to use more than one property at a time. For example, Li et al. (2019) and Wu et al. (2019) used the number of nodes, number of edges, density, degree and modularity together to signal a decrease in bacterial community complexity in response to increasing pollution. Wu et al. (2019) showed that the networks of communities at the more polluted sites in the Jinchuan River (China) had the smallest edges to node ratio, cluster coefficient, average degree and density, which highlighted the reduction in complexity relative to the least and moderate polluted sites. Similarly, in Lake Taihu (China) networks from the least nutrient enriched system contained the highest number of nodes and edges. (Li et al. 2019). However, as in the case of Di Battista et al. (2020), Li et al. (2019) these network attributes did not respond in a dose dependent manner, with both studies observing a lower number of nodes and edges in moderately disturbed systems.

Microbial community networks are dynamic structures that can change over a period of time or along a space gradient, so comparing them can provide insight into how environmental conditions and population variation influence the structure of the microbial population (Pellissier et al. 2018). However, there are some major challenges when comparing multiple systems samples due to being a snapshot in the time and space of the community and because of the absence of a delineated start or end of the ecosystem, they are a continuum. The sample size has to be large enough to capture the entire community of interest, otherwise, the network would not represent all possible interactions and would be misleading (Faust and Shorter 1981; Tylianakis and Morris 2017). It is crucial to control for network properties co-variation during the normalisation step, to enable linear regression measurement between network properties and, to study the residual variation (Pellissier et al. 2018). Although comparing networks usually includes only comparing two relations mapped on the same population during the same period, Faust and Shorter (1981) describe a general way to compare networks using network properties when those networks differ widely in size, type of relation, species of the units, and time and space of the observations (Faust and Shorter 1981). For example, Williams et al. (2014) found that co-occurrence relationships occurring in different environments could be ecologically relevant, such as the relationship between Solirubrobacterales and Acidomicrobiales and the related families, Acidimicrobiaceae and Conexibacteraceae, that are found in soils and the human body (Williams et al. 2014).

Network analysis have an important role in risk assessment, as they allow for modelling outcomes of changes to community structure. A microbial community network is robust when its structure remains stable following a disturbance (Bissett et al. 2013). To test the robustness of a network, a node can be removed or altered to affect the network structure (Barnes et al. 2016). Generally, the removal of the most connected species, hubs, causes more secondary extinctions than a random species (Sheykhali et al. 2019). However, when the system returns to the original state after a perturbation the network is considered resilient (Mandakovic et al. 2018). For example, most denitrifiers are resistant to hydrocarbon pollution because they are insensitive to many toxicants and can adapt to various abiotic factors (Powell et al. 2006).

Co-occurrence networks have great potential as a tool for biomonitoring of aquatic systems, with network properties being potential metrics of change in the underlying community. However, there is still a lack of examples of how to apply this predictive model experimentally in aquatic ecosystems. Most applications of co-occurrence networks have been field-based, where it is often difficult to separate the effect of contamination from other environmental factors, and many studies have contradictory information (Weiss et al. 2016; Hirano and Kazuhiro 2019). Consequently, there is the need for more experimental studies in which confounding factors can be controlled, and the sensitivity of communities and networks to specific stressors can be tested.

The advent of next-generation sequencing has provided an unprecedented capacity to describe microbial communities, but with the limitation that we can now identify many taxa about which we know little of their ecological roles and functions. This provides a challenge to the interpretation of co-occurrence networks because, although we infer interactions between taxa from the networks, there is little empirical evidence to inform (and confirm) the specific nature of those interactions. Further knowledge of microbial interactions under real-world conditions would enhance our ability to understand the particular network elements.

To date, there are only a few examples of co-occurrence network used to test if they respond to a stressor (Pauvert et al. 2019), but none in aquatic environments. Moreover, sample size affects co-occurrence networks (Kolaczyk et al. 2015), a range between 25 and 300 samples per network has been recommended for obtaining reliable networks (Berry and Widder 2014; Faust et al. 2015). Decreasing costs and increasing analytical power of next generation sequencing means that such levels of replication are now feasible, but this remains a potentially limiting factor because of the time and resources required to observe, collect and process the numerous samples necessary, which is well beyond that used in most current monitoring programs (Bohan et al. 2017).

Although in the last decade there has been an increase in the development of technologies, laboratory procedures and bioinformatics pipelines, the workflow should be optimised to reduce handling steps such as extraction, amplification, or sequencing bias (McGee et al. 2019). Moreover, the quality of taxonomy assigned to current OTUs or ASVs references databases varies greatly, some gene databases are well populated and have robust information, and other gene coverage is still incomplete or has incorrectly assigned sequences. Also, it is important to include unknown taxa, as seen by Zamkovaya et al. (2021) the frequent dominance of unknown taxa as hubs stresses the need for further exploration and functional characterisation of all taxa (Zamkovaya et al. 2021) and capture the full taxonomic diversity or the ecological processes being investigated (Zinger et al. 2019). However, co-occurrence network analysis can avoid these issues with taxonomy enabling the possibility to look at the system as a whole and not focusing only on specific nodes.

Networks and their associated properties are new types of indicators that can evaluate the quality of aquatic ecosystems and are promising environmental biomonitoring tools. Co-occurrence network analysis provides a suite of relevant measures that are frequently used in network studies, albeit in various ways, they all use more than one property at a time to understand changes in the network topology.

Two approaches have been found to be most popular. The first approach described the whole network and reveals patterns allowing comparisons between the entire networks. When looking at the whole network topology, the most used properties are the number of nodes, number of edges and the ratio between positive and negative associations that give an idea of the network dimensions. The number of components is used when comparing networks visually giving a piece of information on the integrity of the network, for example, an increased number of components is typical of fragmentation. However, characteristic path length, cluster coefficient and network density are usually used to describe community cohesion and complexity. The second approach used centrality measures such as degree centrality, betweenness centrality, closeness centrality and transitivity, to look at the role of nodes in relation to the rest of the network and allow for comparison of specific nodes, taxa, roles and functions. Moreover, centrality measures are very informative when studying a network’s robustness while removing nodes to affect the network structure. Although there are some examples of network analysis in aquatic ecosystems, they are all field based, therefore, there is still a need to validate the co-occurrence network analysis approach. Consequently, there is an urgency for more experimental studies that are designed to prove whether co-occurrence networks are indeed sensitive to particular stressors or conditions. Owing to the variety of uses, applications and outcomes obtained from previous studies, there is a further need to develop a common standardised approach with clear indication on the sample size, the number of replicates for the same environment and which program to use.

Ultimately, this field is advancing rapidly and we foresee increasing developments in the technical aspects of network analysis with experimental validation to clarify the statistical results, increase our knowledge of the aquatic microbiome and the use of network analysis as a biomonitoring tool.

Data sharing is not applicable as no new data were generated or analysed during this study.

The authors declare that they have no conflicts of interest.

Annachiara Codello received a PhD scholarship from Macquarie University. This research did not receive any specific funding.