Perspectives Paper

From diversity to complexity: Microbial networks in soils

3.1. Preparation of the dataset

In most network construction approaches, associations between microorganisms are identified either from probabilities of co-presence (co-absence) or from pairwise correlations of species abundance in a number of environmental samples. Abundances of microbial taxa in natural microbial communities are usually assessed by amplicon sequencing of highly preserved, phylogenetically informative marker genes, such as the 16S ribosomal RNA for bacteria and the internal transcribed spacer (ITS) for fungi. The sequenced reads are grouped by sequence similarity into operational taxonomic units (OTUs) or recovered from sequenced reads as amplicon sequence variants (ASVs) (Callahan et al., 2017), which can then be classified into taxonomic categories.

Two main characteristics of amplicon datasets are: (i) they are compositional, which means that abundances of taxa can only be interpreted in relative terms (Gloor et al., 2016; Morton et al., 2019), and, especially for soil microbial communities, (ii) taxa are sparsely distributed, which means that datasets contain a high amount ASVs or OTUs that only occur in a fraction of the samples (Alteio et al., 2021). Both of these characteristics pose a major challenge to network construction approaches, as correlations carried out on relative abundance and sparsely distributed datasets can lead to spurious results (Morton et al., 2019; Alteio et al., 2021). In the next subsections we discuss how to reduce the potential bias introduced by these features of the dataset, by applying appropriate data filtering and transforming steps before the actual network construction (Fig. 2).

1. Download : Download high-res image (989KB)
2. Download : Download full-size image

Fig. 2. Workflow for the preparation of a dataset which precedes the network construction. It involves filtering out some ASVs and a data transformation, such as scaling the number of reads a_ij by the total read sum in a sample S_i, or taking a centered log ratio (clr). The transformed dataset contains relative abundance values constrained by a constant (cte) sum T_i, or even inferred absolute abundance values for each ASV in each individual sample.

3.2. Establishing associations

3.2.1. Choice of metric to measure associations and choice of an edge threshold

The most common methods of network construction rely on establishing edges between taxa with significant pairwise associations. In this case the network depends strongly on two key decisions: (1) the choice of a metric to quantify the strength of association's signal; and (2) the choice of a significance threshold - referred to as edge threshold (see BOX1), which serves as a criterium to differentiate a true signal from noise. The strength of associations (1) is most often measured using parametric Pearson or non-parametric Spearman or Kendall-τ correlation coefficients. The latter two are more convenient since they can be used on non-normally distributed datasets, as is the case of microbiomes. However these are not the only possibilities for (1) and measures such as Bray Curtis, Kullback Leibler dissimilarities or even mutual/maximal information can also be used to quantify associations. After the metric is chosen the next step is (2) — the choice of an edge threshold, τ. Before this threshold is applied, the networks are fully connected, with each node (taxa) connected to all other nodes. In these networks the weight of each edge is given by the corresponding pairwise correlation (association) value ρ. To illustrate this, we provide the histogram of Spearman correlation coefficients for our dataset of the upper soil core, which shows that the weights of most of such edges are close to zero and cannot be interpreted as a strong signal of an association (Fig. 3a). The solution therefore, is to eliminate edges with weights below a given edge significance threshold τ. In this way, the threshold is related to the sparsity of the network, where lower values result in very densely connected networks, while high values in networks with fewer edges, see (Fig. 3b). The choice of the threshold value can be either totally arbitrary using a single value for all edges, or it can be evaluated for each edge by comparisons with null models. In the second case, the null model is obtained by shuffling read values among taxa and samples. In other words, a new data set is obtained by randomly selecting values from the original one. The shuffling can be followed by a possible rescaling, due to new total read values in each sample. This procedure is done to break any non-random pairwise association patterns while keeping some of the original characteristics of the data (e.g. sparsity, compositionality and marginal count distributions of taxa). Each pairwise correlation value in the original data set is then ranked against the correlation value of the same pair in the shuffled dataset (Faust and Raes, 2016). This allows one to check the probability of observing such correlation by mere chance (i.e. in similar random data). The p-value determines if the edge should be kept or eliminated. (Fig. 3a) marks with light blue the pairwise correlation values which are left after such a procedure (using 1000 shuffled datasets), and demonstrates that the procedure also works as a sort of edge threshold. Another popular method for edge threshold selection is using ideas from random matrix theory (RMT, mathematical research area that studies matrices with random elements) (Deng et al., 2012). The approach determines at which threshold the generated adjacency matrix no longer has properties compatible to that of a fully random network and uses that as edge threshold. Although the procedure takes away the arbitrariness of the choice, it also selects a hard threshold, and relies on the assumption that biological networks have a completely random structure at the low τ limit.

1. Download : Download high-res image (307KB)
2. Download : Download full-size image

Fig. 3. The choice of the threshold, to establish an association edge, changes the sparsity of the constructed networks. (a) Frequency, N_ρ of the Spearman's rank correlation coefficients ρ obtained for all ASV pairs in the upper soil core (in gray). The selected ρ values by the use of different edge thresholds τ are marked with different colors (red and dark blue). In light blue we show ρ selected as significant by comparison with a null model (for that we use 1000 shuffled versions of our dataset, and p-value = 0.05). (b) Number of positive edges (ρ > 0) in networks constructed with different edge thresholds. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

We illustrate the networks obtained by the use of two thresholds τ = 0.7 and 0.8 in (Fig. 4a, c, e) and (Fig. 4b, d, f) respectively. For both thresholds, networks from upper and lower soil layers, as well as from both soil layers combined present two highly interconnected sets of nodes (modules or clusters): While the two sets are completely disconnected for τ = 0.8, we can detect candidate taxa that connect the two potential communities for the lower threshold. The choice of threshold also changes important network properties, such as mean degree, mean distance and modularity. This demonstrates that the choice of edge thresholds can influence the interpretation of the network structure in terms of its connectivity and the relevance of individual taxa in the system.

1. Download : Download high-res image (2MB)
2. Download : Download full-size image

Fig. 4. Association networks of microbial taxa in a Beech forest soil at different soil depths constructed from Spearman's rank correlation coefficient at different edge thresholds. Shown are networks constructed from upper and lower soil layers combined (a, b), and from each soil layer separately (c–f), for edge thresholds τ = 0.7 (a, c, e) and 0.8 (b, d, f). Notice how the selected value for the edge threshold changes the number of edges in the network as well as the network structure (for details on preparation of the data set see Supplemental material). The network of the combined dataset (a, b) shows not the total 595, but only 461 ASVs shared with either one of the two other datasets (lower panels). The network with the phylogenetic classifications of nodes is also provided in the Supplement. Network visualization was done with graph-tool (Peixoto, 2014).

The advantage of using the strategies described so far is that they involve simple computations. However, as described before they suffer from several limitations: they do not take into account the compositional nature of the dataset and they cannot clearly separate association signals from noise. We can see that the choice of a threshold constitutes a compromise between the detection of spurious links (false positives) and neglecting existent associations (false negatives). The problem of threshold selection is not unique to species co-occurrence networks, but also exist in other areas from inference of gene co-expression to social networks (Couto et al., 2017). The disadvantage of any hard threshold, independent of the selection criteria, is that the thresholding procedure itself may distort the underlying data, strongly affecting the observed network patterns (Cantwell et al., 2020). In particular, as we explain in the next section, this occurs because there is no possible threshold that can distinguish direct from indirect associations. For microbial networks what we often see is that thresholding leads to very dense network structures (‘hairballs’ (Faust, 2021)). All these reasons make it difficult to derive an ecological interpretation from the pattern in the constructed networks. Several network construction tools were developed to tackle the listed issues, we describe some of them in the next section.

3.2.2. Statistical methods that improve reconstruction of ecological networks

In this subsection we mention some methods developed to improve the reconstruction of ecological relationships. However all of these methods depend on the quality of the data and one should not expect signal of interactions in the data if the scale of sampling is inadequate (see Sec. 5.3). In the last decade several ready-to-use tools were made available for network construction. Here we highlight: CoNet, which combines several association measures described in Sec. 3.2.1 (Faust and Raes, 2016); and MENA, which uses random matrix theory to decide on the threshold choice (Zhou et al.; Deng et al., 2012); both are widely used for network construction for microbial soil communities. Another tool which stands out is SparCC due to its ability to treat possible compositionality bias in the data (Friedman and Alm, 2012). It uses log-ratio variance (vlr) to estimate associations, a transformation that presents advantages as it is subcompositionally coherent (see definition in (Table 1). More precisely this transformation allows a consistent analysis of any subset of the data, without any contradictions in the results. The SparCC algorithm works under the assumption of few associations between taxa. The final adjacency matrix is evaluated again either by establishing an edge threshold or a comparison with null models.

To improve our ecological interpretation, it is also of interest to distinguish direct from indirect dependencies between taxa. For example, it is well known that two taxa may show a significant correlation because both are under the influence of a third taxon, although they do not directly interact with each other. The difficulty is that no choice of edge threshold for correlations is capable of distinguishing between direct and indirect associations. To solve this issue a new set of algorithms for network construction was proposed based on probabilistic graphical models (Kurtz et al., 2015; Fang et al., 2017; Yang et al., 2017; Yoon et al., 2019; Jiang et al., 2021). These models represent complex probabilistic relations in the form of diagrams (graphs), where nodes represent taxa and edges represent conditional dependencies between them. Such a graph can be easily constructed for Gaussian multivariate data and for this case it is known as Gaussain graphical model (other names such as partial correlation network or concentration graph are also used). The conditional dependencies are evaluated from non-zero elements of the inverse covariance matrix (Bishop, 2007). It is important to note that strictly speaking for non-Gaussian settings, as is the case for microbial data sets, the relationship between the inverse covariance and conditional dependencies is not known. However, the approach is still used as it is expected to be informative beyond its standard domain (Liu et al., 2009). The inverse covariance is computed with a popular statistical method known as glasso, lasso regularized maximum likelihood estimator. In this case, the sparsity of the network is tuned with a regularization parameter, λ, with small values of λ corresponding to densely connected networks and large values to fully disconnected networks. The method in general selects a range of λ values and generates several networks. The next step is to select the best model, considering its complexity and data fitting ability. The most well known model selection algorithms are BIC, EBIC, Stars. SPIEC-EASI (Kurtz et al., 2015) was one of the first tool packages based on graphical models and therefore is the most well-known. However it is not the only package based on a graphical model, there is a series of other recently developed tools among them: gCoda (Fang et al., 2017), SPRING (Yoon et al., 2019) and BC-GLASSO (Jiang et al., 2021). We also illustrate the outcome of network construction using SPARCC and a SPIECEASI-like (clr transform, glasso and Stars) algorithm in (Fig. 5). For simplicity, we do not show nodes which have zero edges in this image. Again, as in (Fig. 4) it is possible to identify the formation of two clusters, one with the taxa unique to the respective core depth (nodes shown in blue and green), and another with taxa shared by the two cores (nodes shown in red). Despite the differences between methods, they capture similar patterns, although only half of the edges are shared among them, see (Fig. 6). Nevertheless, the ecological processes behind the observed structure remain unclear. Note that these methods perform rather poorly when used on artificial data sets produced even from simple models of population dynamics such as gLV (generalized Lotka-Voltera) when different environmental conditions are used (Berry and Widder, 2014; Hirano and Takemoto, 2019). One of the reasons is that indirect dependencies between taxa can also appear due to environmental factors, for example due to niche overlap. It is possible to introduce environmental parameters as additional nodes in the network (Faust, 2021). Another form to deal with environmental confounders is to use a model based approach as done by mLDM (Yang et al., 2017), a tool introduced recently based on an hierarchical Bayesian statistical framework.

1. Download : Download high-res image (1MB)
2. Download : Download full-size image

Fig. 5. Networks of positive associations for microorganisms in a Beech forest soil at different soil depths constructed by different network inference methods. Networks constructed from upper and lower soil layers combined are shown in (a, b), and from each soil layer separately in (c–f). Networks were constructed by (left panel) SparCC or by (right panel) SPIEC-EASI-like approach (transformed with clr, and then using GLASSO combined with Stars). To improve visibility nodes with no edges are not shown.

1. Download : Download high-res image (324KB)
2. Download : Download full-size image

Fig. 6. The edges detected by different network construction methods, using Spearman correlation (edge threshold 0.7), SPARCC (threshold 0.65) and GLASSO (STARS for model selection) with and the overlap among them.

3.3. Summary of test dataset analysis

We used microbial co-occurrence data from upper and lower soil layers of a European Beech forest to test whether the choices for a specific network inference method (i.e. correlation-based methods with different edge thresholds, (Fig. 4), a method using log-ratio transformations to address a possible compositionality bias, SPARCC, and a method considering indirect associations, GLASSO and Stars; (Fig. 5), or the inclusion of an environmental factor (soil depth) would influence the structure of the constructed network, and the conclusions that can be drawn from it. The datasets of upper and lower soil layers together consisted of 466 ASVs, from which 162 were exclusively occurring only in the upper soil samples, and 183 only in the lower soil sample, and 121 were shared between the two (Fig. S5). We used this taxon-specific information to investigate if and to what extent the resulting network structure was influenced by taxon preference for upper or lower soil, two habitats which differed in many environmental conditions, like temperature, moisture, or substrate input. To visualize this, we categorized the nodes into those which exclusively occur in one of the soil layers after filtering ASVs with low read numbers (Fig. S5), which could be interpreted as ‘specialists’ and those which occurred in both soil layers, which could be interpreted as ‘generalists’. Irrespective of the method used, networks which were constructed including the environmental gradient, i.e. based on the combined dataset from both soil layers, always exhibited two highly connected sets of nodes (clusters), one containing all taxa dominant in the upper and the other those from the lower soil layer. This demonstrates that taxa which co-occur due to similar environmental preferences strongly dominate the network structure, potentially obscuring other ecological relationships.

Ruling out this environmental factor, i.e. constructing networks for each layer separately, revealed an interesting pattern. Here again we see in all methods that two clusters (modules) formed, of which one was dominated by the specialists and the other one dominated by the generalists. This is an interesting finding, and the next step would be to think about its interpretation, or to formulate a testable hypothesis for future research based on this result. Both, however, go beyond the scope of this paper. Here, we just would like to point out that this pattern was conserved, in principle, in all the methods (Fig. 4, Fig. 5). The only observable limitation is that in the most sparsely connected network, i.e. the correlation network with a high edge threshold in the lower soil layer (Fig. 4f), this finding may not be visible anymore owing to the small number of connections. While the coarse pattern seems to be robust and independent of the method, however, we also observed that the majority of links were different between the different methods (Fig. 6). For a more detailed analysis, where pairwise associations between individual taxa are of interest, one should thus consider that methodological choices could influence conclusions at that level. For a display of networks with nodes coloured according to their phylogeny see (Fig. S7 and S8). Interestingly, however, both correlation (τ = 0.7) and SPIEC-EASI like methods identify the same central taxa that intermediated the connection between the modules of the upper soil layer (ASV_32a_pom, a Proteobacteria Rhizobiales). From this small test we conclude that strong patterns in the dataset will be captured by most methods, while the difference in the detailed connection pattern may be quite different, and conclusions on individual connections should thus be done with care.

One interesting question here would be if it was possible to reveal these clustering patterns also using traditional ordination approaches. In principle, ordination-based analyses, such as nonmetric multidimensional scaling (NDMS) or correspondence analysis (CA) can explore and visualize community dissimilarities between samples (i.e. ‘beta-diversity’) or, if turned around, dissimilarities between taxa in their occurrence pattern across samples. One could argue that the latter approach may reveal a potential clustering of taxa in a similar way as co-occurrence network analysis. The reduction of the high dimensional nature of the dataset realized by ordination methods, however, may discard valuable information, affecting the conclusions. Network models, on the contrary, work in the original multidimensional space and therefore are more appropriate to detect such patterns. To test this, we applied taxa-based NMDS to the same dataset we used for the network analysis. While NMDS was also able to separate lower and upper soil specialists in the combined dataset of both layers, the method could not detect a generalist and a specialist-dominated cluster in the dataset of each separate layer as observed by network analysis (see Fig. S9).

4.1. Properties of the nodes

The importance of a species in an ecosystem can be classified according to how many other species it can directly or indirectly influence. In a network of interactions the direct influence of a node is given by its degree (number of edges), see Fig. 7. Node with large degree, or high connectivity, are known as hubs and represent species generalists or ‘key taxa’. While nodes at the network border, classified as peripherals, represent specialists (Olesen et al., 2007). This definition can be extended to encompass more complex network structures and the nodes can be classified in four simple categories: local or global-hubs, connectors or peripherals (Guimerà and Nunes Amaral, 2005; Olesen et al., 2007).

Although this concept to interpret nodes with a high connectivity as ‘key taxa’ stems from plant-animal interaction networks, it was later adopted in a large and still growing number of soil ecological studies to identify ‘keystone taxa’ in soil microbial datasets based on co-occurrence networks (e.g. (Zhou et al., 2011; Lu et al., 2013; Lupatini et al., 2014; Jiang et al., 2016; Ling et al., 2016; Shi et al., 2016; Chen et al., 2018; Qi et al., 2019; Varsadiya et al., 2021; Yuan et al., 2021)). Other attempts to classify keystone taxa in soil microbial co-occurrence networks select nodes with high degree and high closeness centrality, but low betweenness centrality (Banerjee et al., 2018, 2019; Lin et al., 2019), guided by a theoretical analysis of co-occurrence networks constructed from generalized Lotka-Volterra equations (Berry and Widder, 2014). It must be noted that one of the selection criteria in this approach, low betweenness centrality, differs fundamentally from the selection criteria of high connectivity defined in (Olesen et al., 2007), so that the two approaches would likely lead to the identification of different keynote species when applied to the same dataset.

Most important, however, both (and other) approaches to identify ‘keystone species’ from soil microbial datasets currently suffer from a lack of conceptual justification of why a node with a certain topology, for example a highly connected (‘hub’) node, should indicate a ‘keystone’ role of the respective taxa in soil microbial co-occurrence networks (Röttjers and Faust, 2019). A ‘keystone species’ is defined as ‘a species whose impact on its community or ecosystem is large, and disproportionately large relative to its abundance’ (Power et al., 1996). This means that the removal of such a species will have a detrimental effect on the community or ecosystem, a behavior which has been computationally shown when removing hub species in ecological networks (Dunne et al., 2002; Memmott et al., 2004). A requirement for transferring this conclusion from the constructed network behavior to ecosystem behavior, however, is that the edges in the constructed network represent vital and essential links between species in nature. While an edge in an interaction network, such as the plant-pollinator networks described in Olesen et al. (2007), reflects with a certain likelihood an essential interaction between two species, so that the removal of a hub species may indeed lead to the loss of other species, this cannot be assumed for soil microbial co-occurrence networks. Even if the simulated removal of a node in a microbial co-occurrence network affects the computed network structure, it would not need to have the same effect in nature, as species who co-occur do not necessarily depend on each other.

4.2. Average path length of the network and degree distribution

The architecture of biological networks is often far from random. In general, interaction networks, such as food webs and mutualistic networks are known to have broad degree distributions, with few very connected nodes (hubs) and the majority of species that interact only with few others (Camacho et al., 2002; Montoya et al., 2006). This degree distribution also leads to short path length between nodes, meaning that randomly chosen species have only few edges between them (Williams et al., 2002). The questions are how and why such types of architecture develop for microbial communities and their relation to stability or community assembly mechanisms. In the case of co-occurrence networks there is need to disentangle possible reasons behind the observed patterns, such as evolutionary history, dispersal dynamics or negative/positive biotic interactions (Weiher et al., 2011; Morueta-Holme et al., 2016; Goberna et al., 2019). Relevant questions that still have to be approached are for example: how relative distances between taxa in the network are related to either phylogenetic similarity or even to physical distance if environmental gradient is present.

4.3. Cluster (network community) detection

Ecological networks can be composed of hundreds, or sometimes even thousands, of nodes. It is natural, therefore, to wish for identification of representative high-level subsets or groups to organise this complex information. However, we do not want any partition, but the one that can inform us of underlying functional, evolutionary or environmental units within these ecosystems. This is an easy procedure if we have metadata on such groups beforehand, e.g. division of microorganisms in groups predominant in certain environments (Eldridge et al., 2020). However, in most cases we lack this information and the underlying partition is precisely what we would like to learn about the system. The problem therefore is to determine node groups based on network topology, for instance by checking the density of connections among them (Newman, 2010). This is known as cluster detection or community detection. Note that this term is used for ‘network communities’, which consist of several similar nodes, not to be confused with microbial communities. This task consists of: (1) finding a natural division in the network, without being able to look at its structure; and more importantly (2) evaluating if these patterns are statistically/ecologically meaningful or if they can appear by mere chance.

Certainly the most popular algorithm developed for cluster detection (1) relies on modularity maximization, which formulates the task of community detection as an optimization problem. This is done by using a mathematical function known as modularity to assign a score to each division of the network (Newman and Girvan, 2004; Newman, 2004, 2006). Although the method has attracted substantial popularity due to its simplicity, it has drawbacks which are often ignored in applications. While the method works well when a clear division exists, it is also known to find spurious partitions when there are none to be found, e.g. in completely random networks (Guimerà et al., 2004). This conflicts with our objective (2) of selecting only statistically meaningful partitions. Moreover, a comparison with a null model cannot solve issue (2) since it can only rule out or accept the null model (random network structure), but cannot assess if the found partition is in fact statistically significant. This is because finding enough evidence to reject a null model is a much weaker requirement than having enough evidence to accept a particular network partition.¹ Due to these problems several other methods have been proposed to extract meaningful communities from networks (Fortunato, 2010; Fortunato and Hric, 2016). We highlight the inference methods based on stochastic block models (SBMs) (Karrer and Newman, 2011; Abbe, 2018; Peixoto, 2019). Their main advantages is their ability to distinguish signal from noise, taking care of (2) in a principled manner, and possess also the ability to detect not only patterns based on assortative structures (Newman, 2004), but others, such as core periphery, or self-similar hierarchical patterns (Fig. 7).

The kinds of patterns identified by cluster detection have been proved useful for understanding ecological networks: Structural assortative as well as nested patterns have been linked to stability in interaction networks (Thébault and Fontaine, 2010; Grilli et al., 2016). In association networks clusters have mainly been used to identify and separate different niches occupied by microorganisms.

4.4. Network stability

Over the past decades an increasing number of studies have approached the question of ecosystem stability by analyzing the architecture of biological interaction networks (Montoya et al., 2006; Okuyama and Holland, 2008; Bascompte, 2009; Landi et al., 2018). The stability is characterized by how systems respond to disturbances, and whether or not they can return back to their original function after being perturbed. It can be evaluated with the use of different criteria, measuring the recovery process (resilience), endurance (robustness) or the flexibility (persistence) of an ecological system (Hodgson et al., 2015; Landi et al., 2018). Although sometimes these terms are used interchangeably, what is important to note is that the type of perturbation and the type of analysis may depend on the research question. For example one can investigate how a community endures species extinctions (robustness). Alternatively, one can observe how quickly the multi-species dynamics recovers to a decrease of species abundances or to invasion of a new taxa (resilience).

A significant fraction of the studies on stability are theoretical and use numerically generated datasets (e.g. gLV, generalized Lotka-Voltera model). They address a long standing debate in ecology on the relationship between the ecosystem's complexity (e.g. number of connections in an ecological network) and its stability. The controversy in the field started around the 1970s when theory showed a conflict between the two aspects (May, 1972, 1974; Sales-Pardo, 2017). Currently we know the network properties that foster stability, but only for some types of networks. For example, nestedness and connectivity seem to promote stability in mutualistic (Thébault and Fontaine, 2010) and resource competition networks (Wei et al., 2015), while for food-webs a compartmentalized structure, low connectance seem to bring stability (Neutel et al., 2007; Gross et al., 2009). In the presence of both positive and negative interaction evaluating the stability is more difficult. It is speculated that the presence of mutualistic interactions may destabilize ecosystems, since they generate co-dependencies (Coyte et al., 2015). Trophic interactions, on the contrary, were shown to promote stability due to negative feedbacks (Coyte et al., 2015). Due to the difficulties presented above, claims about stability are even more unreliable for microbial co-occurrence networks, where there is not only the possibility of different interaction types but the role of links is less clear.

On the experimental side, studies have focused on characterizing how the ecological networks and ecosystem function change in reaction to a perturbation, or along gradients of environmental conditions. Some of these studies showed that connectivity, density of the links and transitivity decrease in disturbed environments (Karimi et al., 2016, 2017; de Vries et al., 2018). Going one step further, some studies used soils that were exposed to different treatments or perturbations to investigate a possible link between the network structure and the functionality of microbial communities, or even ecosystem multifunctionality (Wagg et al., 2019; Yuan et al., 2021; Qiu et al., 2021; Jiao et al., 2022). However these studies based their conclusions primarily on correlations between network complexity and functionality measures, and were thereby unable to provide any mechanistic foundations for their claims. Again, we recommend to be cautious to draw causal relationships from correlative pattern, which we think is especially problematic for correlations between two highly aggregated measures, such as network complexity and soil or ecosystem multifunctionality, as they naturally may share a high number of co-correlating factors and interdependencies. Therefore, instead of trying to correlate given functions with network patterns, we suggest to use studies encompassing different treatments or perturbations to test hypotheses that for example link biodiversity and complexity. Such hypotheses can be formulated based for example by randomly (or systematically) removing nodes from an unperturbed network and recording how the connectance of the network changes. These hypothesis can then be validated empirically by comparing it to networks constructed from observational data of a different treatment (where a perturbation or stress is applied), checking if the loss of diversity affects the network in a similar way. Such comparisons can also allow to analyze if perturbations or stress leads to loss of edges without a loss of corresponding taxa (Valiente-Banuet et al., 2015) or there is a formation of new patterns. Formulation of hypotheses on the resilience of a system can be an important way to validate network models inferred by the current methods. For example, tracking the presence of known keystone species and how the network structure is affected by them can be of interest. Finally, we would like to mention the appeal of using co-occurrence networks as indicators of ecosystem's quality and function as proposed by Karimi et al. (2017).

5.1. Specific challenges in soil microbial ecology

Although there are several possible drivers behind the formation of microbial communities, species-species interactions have taken the spotlight of current research on taxa co-occurrence patterns. Despite its appeal, the detection of interactions from this type of data has proven to be extremely challenging. Moreover, it is still unclear whether our limitation in spotting species-species interactions can be attributed to an inefficiency of current inference methods, or to the weakness of interaction signal in abundance (co-occurrence) data (Blanchet et al., 2020).

Clearly, besides the potential signal of species interactions, association networks contain other valuable information about relationships within an ecosystem and can help us to understand ecological processes behind the high microbial diversity characteristic to soil. However it is important to emphasize that any network model and its analysis is only as good as our interpretation of it. Therefore, it is desirable to include an interpretation of all associations represented as edges to the most possible extent. This is particularly hard for soils, since they are characterized by their immense biological diversity and heterogeneous physical structure. Here we summarize the important soil characteristics that may impact network construction, interpretation and analysis, and subsequently discuss how we can use experimental design and statistical methods to improve our ability to reconstruct ecological relationships.

5.2. Soil microbial networks: what type of associations should we expect?

Soil provides a unique environment for microbes, because it is spatially structured and chemically and physically heterogeneous across several scales. Within the complex soil matrix, numerous edaphic properties impact the structure of the microbial communities therein, including pH, salinity, temperature and moisture (Frindte et al., 2019; Zheng et al., 2019). Additionally, temporal dynamics such as weather patterns, root exudation and other seasonal inputs of organic material may influence the structure and activity of a microbial community (Kuzyakov and Blagodatskaya, 2015; Chernov and Zhelezova, 2020). As a result of soil heterogeneity and temporal fluctuations, it is critical to consider that sequence-based approaches reveal only a snapshot of the microbial community present at a given time. Together, these properties of soil as a habitat for microbes influence the possible scope for interpretation of microbial association networks.

Heterogeneity of soil samples is a major source of concern in the development of reliable inference methods to uncover ecological relationships. We see the heterogeneity impacting the sampling and inference in two different ways. First, it can lead to hidden physical or chemical differences between apparently similar samples (Armitage and Jones, 2019). Since environmental drivers lead to strong differences in the community structure, such drivers may dominate the detected association signal. In this case, additional environmental measurements can help to interpret the observed patterns. The second concerns the large volume of each sample used in the analysis. Nucleic acids extracted for sequencing studies are derived from amounts of soil ranging from 250 mg to 2 g of material. Very often, they are aliquots from an even larger soil volume, for example several soil cores that have been mixed and homogenized to be representative for a certain field site or plot. Soil microbes, however, interact with each other at the scale of tens of μm (Raynaud and Nunan, 2014). It is obvious that co-occurrence of microbial taxa across soil samples taken at a scale which is by orders of magnitude larger than the scale at which they interact, does not necessarily mean that they interact with each other. In line with this, community assembly theory suggests that the larger the sampling units scale the more environmental filtering effects dominate species co-occurrences (Kraft et al., 2015). Although the mix of ecological factors operating at larger sampling scales obscures causal distinctions, studies have shown that networks generated with this type of data hold predictive value for central ecosystem functioning (Shi et al., 2016; Wagg et al., 2019). We should, however, lower our expectations to capture a clear signal of microbial interactions, since even small soil samples may contain numerous metacommunities - small consortia linked by dispersion and diffusion (Armitage and Jones, 2019). Two microbial taxa can be separated by unconquerable distances even in a small piece of soil given its microscale physical structure and the spatial heterogeneity of microbial niches. In view of this complexity it is controversial if in such samples biological interactions can leave any signal at all. For example, models suggest that negative interactions, such as competition, can be completely hidden at scales much larger than the scales of the interactions themselves (Araújo and Rozenfeld, 2014). In summary, the associations described in co-occurrence networks from standard soil samples have to be carefully interpreted. Such networks are expected to be dominated by associations driven by the environment, however the obtained network structure may be a result of several unknown factors. In relation to interactions among organisms, especially the signal of competitive interactions will be obscured. In the next section we give suggestions on how to disentangle the different factors behind a given co-occurrence dataset and suggest alternative or complementary experiments.

5.3. Reconstructing networks of ecological relationships

5.3.4. Descriptive, exploratory or hypothesis-based approaches?

Recently we have observed an increase in the number of articles which use microbial networks as a form of data visualization for co-occurrence patterns in soil. These studies define the edges of networks as ‘co-occurrences of microorganisms’ (or even incorrectly as microbial interactions), and have a mainly descriptive character. Without presenting hypotheses of biological mechanisms behind the observed statistical patterns in co-occurrence data, such studies are of little interest to soil microbial researchers. In the previous section we presented several research questions which can be addressed with network analysis. Addressing these questions using the constructed network can be a simple form to extend and improve such descriptive studies. Networks should work as a tool to help us understand the processes that structure microbial communities, and therefore the interpretation of observed patterns and creation of hypotheses from these patterns are not optional and should be included.

Finally, we shortly discuss possible ecological hypotheses that might be tested or used in microbial network analysis. As described in Sec. 4.1, identification of potential keystone species or functional components of the network can help us to forecast structural shifts in the network in response to environmental perturbations, see (Fig. 8c). Another possible research direction is to compare the network architecture shaped by particular conditions, e.g. the available substrate, see (Fig. 8c). In this context, networks can identify taxa with potential similar functionality, since those are expected to compete and be negatively associated across communities (Brown and Wilson, 1956). In the same way we can formulate hypotheses on functionally complementary species, expected to be present in the degradation of complex substrates (Lindemann, 2020). This can be done using networks constructed from datasets extracted from different environmental conditions, expecting a change in predominance of facilitation to competition with increasing nutrient availability (Hoek et al., 2016). These nutrient-based community responses can be combined with a stress response (stress imposed by an environmental factor). What is expected then, is that high stress would reduce investment into competitional traits, while facilitation would enable stress mediation, equivalent to what we know from plant ecology (Bertness and Callaway, 1994; Hammarlund and Harcombe, 2019; Piccardi et al., 2019).

1. Download : Download high-res image (810KB)
2. Download : Download full-size image

Fig. 8. Use of network analysis in exploratory vs hypothesis driven research. (a) Taking into account that community composition is shaped by underlying ecological mechanisms driven by the physical environment or by interactions among taxa. (b) In general, the measured co-occurrence patterns are used to construct association networks, which are different from the interaction networks. (c) Hypothesis-driven research is key to design comparative experiments aiming to isolate and contrast the community assembly mechanism of interest. In particular, established ecological theories can complement hypothesis generation.