A framework for comparing microbial networks reveals core associations

Microbial network construction and analysis is an important tool in microbial ecology. As microbial interactions are challenging to infer experimentally, such networks are often constructed from statistically inferred associations and may not represent ecological interactions. Hence, microbial association networks contain a large number of errors and their derived properties do not necessarily reflect true community structure. Such errors can be identified with the use of appropriate null models. We have developed anuran, a toolbox for investigation of noisy networks with null models, for identification of non-random patterns in groups of association networks. This toolbox compares multiple networks to identify conserved subsets (core association networks, CANs) and other network properties that are shared across all networks. Such groups of networks can be generated from a collection of time series data or from cross-sectional sample sets. We use data from the Global Sponge Project to demonstrate that different orders of sponges have a larger CAN than expected at random.

Introduction 48 across networks not driven by similarity in species composition. The edge intersection of 49 networks can become large even for completely random networks due to similarities in 50 the abundance data. These similarities are relevant when ecosystems consist of a 51 number of "core" genera that are found in most samples of the ecosystem. Such core 52 microbiomes have been identified for the human gut (Falony et al., 2016), the oral 53 microbiome (Zaura et al., 2009) and the coral microbiome (Ainsworth et al., 2015). Yet,54 presence of a core microbiome does not necessarily imply presence of a core microbial 55 interaction network, since the same microorganisms may interact in different ways, 56 fluctuations of their abundances may not be driven by interactions, or the data may be 57 too noisy to infer associations. This raises the question whether interactions are 58 preserved across different realisations of an ecosystem, whether they are preserved 59 within sub-groups of the ecosystem or whether they are unique (Bashan et al., 2016). 60 Since null model analysis can distinguish random from significant intersection sizes, it 61 can test for the presence of one or more core association networks (CANs) within a 62 group of association networks. Core associations can then be further explored to check 63 whether they represent ecological interactions. In this way, null model analysis is a step 64 towards answering whether interactions are universal.

66
For investigation of non-random network properties, several network null models 67 have already been implemented in previous research or in software. Prior work 68 demonstrated that many network properties can be highly correlated, but may represent 69 different underlying aspects in network topology. For example, Valente et al. (2008) 70 found strong correlations between degree and several other centralities, including 71 eigenvector and betweenness centrality in their analysis of 62 sociometric networks 72 (Valente et al., 2008). For bipartite networks, Dormann et al. (2009) demonstrated that 73 the degree and edge density were strongly correlated to species number, while the 74 cluster coefficient and connectance increased as the number of observed edges per 75 species (sampling intensity) increased (Dormann et al., 2009). Consequently, Dormann 76 3/25 et al. (2009) implemented null models in the bipartite R package to correct for sampling 77 intensity. The BiMat MATLAB package uses a similar null model strategy to estimate 78 statistical significance of network properties including nestedness, modularity and 79 module structure (Flores et al., 2016). Connor,Barberán and Clauset (Connor et al.,80 2017) use Chung-Lu and Erdős-Rényi models to identify whether average path length, 81 modularity, diameter and the clustering coefficient were significantly different for 82 observed networks compared to networks generated with these models. In each of these 83 examples, a network property of interest could be interpreted more meaningfully after 84 comparison to random networks. In contrast to studying one network in particular, we 85 adopted a null model strategy to specifically address the existence of conserved In this manuscript, we introduce a null model strategy based on shuffled networks 92 (Fig. 1). We illustrate the power of this strategy by identifying non-random CANs in 93 sponge microbial networks. For these networks, we show that the CANs represent 94 biologically relevant group-specific associations. We also demonstrate how 95 degree-preserving randomized networks help identify betweenness centrality scores that 96 result from the network's connectivity pattern rather than from node degree 97 distribution alone. These strategies have been implemented in a software toolbox that 98 evaluates the significance of network or network property comparisons.

101
Null model toolbox 102 We have developed a software toolbox, anuran, (a toolbox with null models for 103 identification of non-random patterns in association networks) that generates random 104 networks and assesses properties of these networks. Three types of networks can be 105 4/25 generated in the current implementation: completely randomized networks, degree-preserving networks and a variation of both networks that keeps a fraction of the 107 edges fixed. For clarity, the completely randomized networks and degree-preserving 108 networks without fixed edges are referred to as negative control networks, and the same 109 networks with fixed edges are referred to as positive control networks in the remainder 110 of the manuscript. For the completely randomized model, a network is initialized with 111 the same nodes as the input network. Edges are then added randomly until the total 112 edge number is equal to the number of edges in the input network. For the 113 degree-preserving model, edges are swapped rather than removed and added back to the 114 network, so that two edges (a,b) and (c,d) become the new edges (a,c) and (b,d). Hence, 115 the model preserves the degree distribution found in the input network and each node 116 has the same degree as it has in the original network, but other centralities such as the 117 betweenness centrality can change. The user specifies both the number of random 118 networks generated for each network (by default 10) and the number of sets of these 119 networks (by default 50) that are sampled to calculate set sizes. The toolbox can generate random networks with a fraction of fixed edges, which 122 serve as a type of positive control. First, a fraction of edges is extracted from the total 123 union of edges across all networks. For fully randomized networks, these edges are first 124 added, then edges are added until the total number of edges in the original network is 125 reached. For the degree-preserving randomized networks, negative control networks 126 (with preserved degree) are first generated. Then, for each edge in the fixed core, the 127 algorithm attempts to find two edges that can be swapped so the fixed edge is created. 128 If this fails, a random edge is deleted and the fixed edge is introduced, so the degree is 129 not exactly preserved. To swap the edges successfully, it is necessary that each of the 130 nodes participating in a fixed edge has another edge not part of the fixed core. As a 131 result, the degree distribution can change significantly for networks where nodes in the 132 fixed core are disconnected, or where the fixed core is very large compared to the  It is possible to include nodes without significant associations in the network file as 137 5/25 disconnected nodes (orphan nodes). This can be done by simply supplying the network 138 file with the orphan nodes included as nodes without any edges. In this case, the 139 random model more accurately reflects a situation where associations are randomly 140 selected from the number of taxa in the ecosystem. However, the degree-preserving 141 networks will not change, since an orphan node has a degree of 0 and will retain a 142 degree of 0 in the null model. The inclusion of orphan nodes leads to different estimates 143 for set sizes for the random model that may lead to an overestimation of the significance 144 of a CAN. We expect that many nodes present in the abundance file will not be 145 prevalent enough, not abundant enough or not variable enough to ever have a significant 146 association in the first place. Therefore, we ignored the presence of disconnected nodes 147 in our case study.

149
The toolbox has been implemented in Python 3.6 and consists of both an application 150 programming interface (API) and command line interface (CLI) that carries out a 151 complete analysis. Functions from the API can be used to integrate the toolbox into 152 other projects and pipelines and have been structured in a modular manner so 153 customization of the toolbox is straightforward. For example, users can incorporate any 154 centrality or network property implemented in another library, e.g. NetworkX, without 155 making large changes to the toolbox (Hagberg et al., 2008). Currently, the CLI pipeline 156 assesses set sizes, (rank-transformed) betweenness, degree and closeness centrality scores 157 and several network-level properties: degree assortativity, connectivity, diameter, radius 158 and average shortest path length ( Fig. 1). NetworkX implementations of these 159 centrality calculations were used (Hagberg et al., 2008).

161
The software uses a set-of-sets approach to identify CANs. A set is a specific 162 collection of edges, such as the intersection set, which is the collection of edges present 163 across multiple networks. The CANs are identified as differences of specific intersection 164 sets. Hence, the toolbox specifically identifies sets and sets-of-sets that are likely to be 165 of interest for microbial association networks. These sets represent collections of edges 166 that are only present in one specific fraction of networks and distinguish between 167 less-conserved and more-conserved edges. diagram plus the other colours inside this part. To obtain the difference of the 171 intersections, the set that includes one or more additional networks is subtracted from 172 the intersection set that includes fewer networks. These sets are referred to as 173 combinations of intersections with fractions or integers. Therefore, the intersection 0.5 174 refers to all intersections of 50% of the networks, while the intersection 3 refers to 175 intersections that contain edges present in at least 3 networks. Similarly, set-of-sets are 176 identified by a combination of intersection numbers: the set-of-sets 6 → 10 refers to the 177 difference of intersection 6 and intersection 10 and therefore contains no edges present 178 in at least 10 networks. For most analyses, the difference of intersections is preferred 179 over intersections since the intersections are nested. By taking the difference, it is 180 possible to distinguish between more and less conserved associations.

182
The equations for differences and k -intersections for groups of n networks are given 183 below. The equations only refer to edge sets E, so they do not apply to numbers of 184 matching nodes. The difference is the union of all sets D i for 1 up to n networks, where 185 the sets D i contain all edges x present in an edge set E i but not in the union of all 186 other edge sets:  Hence, a k -intersection is the union of all intersections S I for I in P n k : Since edges present in at least k networks but not in m networks represent 198 less-conserved edges, the difference of the intersections is calculated to distinguish 199 between less-conserved and more-conserved edges. The difference of two intersections k 200 and m, with S I and S J defined identically to S I in the equation above, is then given To compare observed set sizes to set sizes of random networks, the Z-score test is 203 carried out, which identifies set sizes in the input networks that are outside the range of 204 indices (from 1 to the network number) and centralities or other network properties.

225
The results of these analyses are exported to tab-delimited files so they can be further 226 analyzed and visualized in the user's preferred statistical environment. rarefied to even depth. After rarefaction, the abundance data were first filtered for 20% 238 taxon prevalence across all samples, then once more to ensure 20% prevalence across the 239 different orders. Counts for removed taxa were retained to preserve the sample sums.

240
After excluding host orders with fewer than 50 samples, 10 orders remained. CoNet was 241 then used to infer association networks (Faust & Raes, 2016

269
We analyzed 10 sponge order-specific networks that we inferred from Sponge  Intersection differences up to 6 networks were significantly larger than differences 289 generated from the randomized and degree-preserving negative control networks 290 (p<0.0001) (Fig. 3). However, the positive controls with a core conserved across 50% of 291 networks had a much larger set size at 5 networks. Therefore, the CAN was constructed 292 from all associations present in 3 out of 10 networks (Fig. 3) In addition to a CAN, we investigated whether the networks contained taxa with 305 consistently lower or higher centrality scores compared to the randomized networks 306 (Additional File 1: Figure S7).  (Burgos et al., 2007). Hence, robust estimates of these properties can yield valuable 322 information on community structure.

324
We chose to use two types of network null models that represent two extremes in 325 terms of constraints. The randomized null model is not constrained in terms of network 326 structure (apart from node and edge number), while the degree-preserving models may 327 be overly constrained especially if the degree centrality is a meaningful representation of 328 a biologically relevant property, such as a taxon's generalist lifestyle. Moreover, these 329 models make no assumptions on the nature of associations. Therefore, it is unknown 330 whether an association in a CAN reflects a biotic interaction, as these associations can 331 also result from similar taxon responses to the environment or to other organisms. The 332 effects of biotic interactions can be better studied through other methods, for instance 333 joint species distribution models Tikhonov et al., 2017).
addition to the effect described in the pseudo-null model. A comparison to these 349 network models therefore addresses whether properties are significantly different 350 compared to networks generated according to specific rules of network growth, but it 351 cannot address the non-randomness of network properties. Since this toolbox has been 352 developed to find non-random trends in association networks, we chose network null 353 models that make no assumptions on network growth (e.g. preferential attachment) 354 because 1) we cannot be sure that such assumptions hold for interaction networks and 355 2) these null models do not take additional processes into account (such as 356 environmental influence) that likely shape association networks.

358
As one of the null models in anuran preserves the degree distribution, comparisons 359 to networks generated from those null models support statements on correlations 360 between degree and other centralities. While we did not fully explore the effects of other 361 properties, such as the fraction of realized edges or network topology, these topics have 362 been discussed previously in methodological studies on ecological and social networks 363 (Dormann et al., 2009;Faust, 1997;Valente et al., 2008;Valente & Foreman, 1998).

364
However, we do want to emphasise that they deserve similar attention in the context of 365 microbial association networks. For example, Agler et al. (2016) defined hub taxa as 366 those taxa that had both higher closeness, betweenness and degree centrality than other 367 taxa (Agler et al., 2016), but such measures may be strongly correlated in dissortative 368 networks, where nodes with high degree are more likely to connect to nodes with low 369 degree (Goh et al., 2003). Our analysis of centrality rankings further supports the 370 observation that betweenness and degree centrality can be correlated, as we found that 371 taxa did not have betweenness centrality rankings significantly different from 372 betweenness centrality rankings observed for the degree-preserving random networks.

374
On networks constructed from sponge order-specific taxon abundances, the 375 set-of-sets approach identified a large CAN. This CAN was significantly different from 376 13/25 the CAN observed for random networks. We suspected based on prior work that this 377 large CAN could arise from a major partition in sponge symbiotic relationships: 378 sponges tend to either have high or low microbial abundance (Gloeckner et al., 2014). 379 Several taxa have been identified as indicators of this divide and many of those 380 indicators were also found in the CAN . Although 381 HMA-LMA status is not strictly phylogenetically conserved across most sponges 382 (Gloeckner et al., 2014), associations between taxa that relate to this status appear to 383 be conserved across at least a subset of sponge orders ( Figure 5). Hence, the set-of-sets 384 analysis suggests that some of the dynamics responsible for the HMA-LMA discrepancy 385 are shared across different orders of sponges. Therefore, the flexible CAN threshold that 386 we adopted has the potential to identify group-specific dynamics that can reveal more 387 about community structure. construct correlations across ordered networks. Our set-of-sets approach uses null 395 models to find conserved patterns across groups of networks. We demonstrate that such 396 a set-of-sets approach, by using an alternative null model strategy, can also detect 397 redundancy between network properties such as degree and betweenness centrality.

398
Therefore, anuran is one of the first dedicated tools for meta-analysis of noisy networks 399 with null models. We expect this null model suite to be a valuable benchmarking tool in 400 the analysis of microbial and other networks.  Figure 1. The anuran pipeline. a Networks are imported by the user. These networks can be ordered, and multiple groups of networks can be imported at the same time. b Random networks are constructed for each of the imported networks. These can be fully randomized, or they can preserve the degree distribution of the original network by swapping edges (highlighted in yello. They may also contain a synthetic core. c The toolbox returns sizes of sets (e.g. intersections of 3 networks) and sets of sets (e.g. the difference between the intersections of 2 networks and the intersections of 3 networks). These sets measure the overlap between specific numbers of networks. By comparing set sizes for observed networks and random networks, the significance of the set sizes can be assessed. d Network-level and node-level properties are computed for both the random networks and the observed networks, which allows assessing the significance of these properties within a group of networks and between multiple groups of networks.

16/25
A B Figure 2. Comparison of FlashWeave and CoNet networks on sponge orderspecific networks. Host order-specific networks were generated with both CoNet and Flashweave from samples collected for the Sponge Microbiome Project. The labels on each figure refer to the sponge order. a Sizes of the FlashWeave and CoNet networks. b Difference and intersection of the FlashWeave and CoNet networks. The different bars represent the observed set size for the observed data (Input), the fully randomized networks (Random) and the degree-preserving randomized networks (Degree).

Negative control
Positive control -prevalence Input Random Degree

Random -10%
Degree -10% Random -50% Degree -50% Figure 3. Set sizes across networks for 10 sponge networks and randomizations of these networks. The set size is the number of edges present in a particular number of networks. The set size is shown for a set-of-sets of network intersections, meaning that the set-of-sets 4 → 6 is calculated as the number of edges in 4 or 5 networks with all edges in at least 6 networks removed. Each network was generated for a different host sponge order for which at least 50 samples were available. These networks were then randomized either with the same degree distribution (Degree) or without this distribution (Random). Moreover, each of the networks was randomized with preservation of a part of the input network for a subset of the randomizations as a positive control. Hence, the positive control Degree networks are randomized versions of each input network with 20% of the union of associations present in at least 20% of observed networks or at least 50% of observed networks. Error bars represent the standard error across different combinations of random networks. For sets of edges present in up to 6 networks, the set size of the input networks deviates significantly from the set size from those of random networks with or without degree preservation.   and were constructed using samples from a single order of sponges. Associations present in at least 3 networks were included in the CAN. Node colour was mapped to phylum and the manta clustering algorithm was used to identify clusters in the network. Figure 5. Confidence intervals of betweenness centrality rankings for three taxa found in sponge order-specific networks. Estimated upper limit and lower limits for the 95% confidence intervals of betweenness centrality computed for each taxon across a group of networks and null models generated from these networks. Taxa that consistently have high centralities have narrower confidence intervals and are therefore located in the upper right of the plot, while taxa with consistent low centralities are located in the bottom right. Null models either had the same degree distribution (Degree) or were fully randomized (Random). The shown taxa were selected because they had a p-value below 0.3 for any of the permutation tests comparing closeness, betweenness or degree centrality rankings. The taxa labelled SAR202, iii1-15 (a clade of Acidobacteria) and Alphaproteobacteria had p-values of 0.119, 0.188 and 0.129 for the degree centrality tests compared to the fully randomized models.