The Self-Organization of Marine Microbial Networks under Evolutionary and Ecological Processes: Observations and Modeling.

Evolutionary and ecological processes are primary drivers of ecological network constrictions. However, the ways that these processes underpin self-organization and modularity in networks are poorly understood. Here, we performed network analyses to explore the evolutionary and ecological effects on global marine microbial co-occurrence networks across multiple network levels, including those of nodes, motifs, modules and whole networks. We found that both direct and indirect species interactions were evolutionarily and ecologically constrained across at least four network levels. Compared to ecological processes, evolutionary processes generally showed stronger long-lasting effects on indirect interactions and dominated the network assembly of particle-associated communities in spatially homogeneous environments. Regarding the large network path distance, the contributions of either processes to species interactions generally decrease and almost disappear when network path distance is larger than six. Accordingly, we developed a novel mathematical model based on scale-free networks by considering the joint effects of evolutionary and ecological processes. We simulated the self-organization of microbial co-occurrence networks and found that long-lasting effects increased network stability via decreasing link gain or loss. Overall, these results revealed that evolutionary and ecological processes played key roles in the self-organization and modularization of microbial co-occurrence networks.

1. Introduction

The alterations of species interactions in ecological network, whether direct and indirect, are giving rise to cascading extinctions [1,2]. Thus, there is an urgent need to explore the determinants of species interactions and the ways that they drive assemblage dynamics of ecological networks. Direct species interactions are widely considered to be constrained by evolutionary and ecological processes [3,4,5,6,7,8]. For examples, phylogeny and optimal niches or functional traits, the proxies of evolutionary and ecological processes, respectively, are key predictors of species interactions [9,10,11,12,13]. However, the question of which process plays the dominant role in constructing species ties is still understudied. Additionally, few studies pay attention to the relative contributions of these processes to indirect species interactions within ecological networks.

Ecological networks are inherently hierarchical [14], and thus the assembly driven by evolutionary and ecological processes can be characterized from various network levels, such as nodes, sub-networks (e.g., motif [15] and module), and the whole network. The metrics at node level, e.g., the degree [16] and within- and among-module connectivity [17], could measure the topological role of species in ecological networks [18,19,20]. However, less is known on how such topological properties are jointly driven by the phylogeny and niches of these species. At the sub-network level, motifs are common network backbones, and their elements may have similar traits or shared evolutionary histories [15,21,22]. Modules are purely identified by mathematical methods, such as the algorithms of greedy modularity optimization [23] and short random walks [24], but constrained by evolutionary and ecological processes revealed by previous reports on pollination networks, seed-dispersal networks and microbial networks [17,20,25,26]. At the whole network level, both processes are observed as primary drivers that shape microbial networks [27,28,29], but how these processes shape network global properties, such as self-organization and modularity, is less reported.

Such self-organization and modularization of ecological networks could be essentially revealed by establishing a mathematical model and simulating network assembly. In order to construct ecological networks, there are the three key features that must be re-emphasized from previous studies [30]: (i) the aggregations of nodes are dynamic; the network can gain or lose them, and they can connect or disconnect with each other [31], meaning that species have dynamic connections with other species. (ii) The nodes harbor specific properties, such as phylogeny and niches. (iii) External forces such as environmental disturbance can also act on network structures. Moreover, we expect that evolutionary and ecological processes should further take into account the self-organization of ecological networks. It is critical for us to understand how ecological networks respond to environmental changes and access network stability beyond posteriori knowledge.

Here, our objective is to study the evolutionary and ecological processes that shape marine microbial co-occurrence networks across multiple network levels. Based on global ocean metagenomic data from the TARA project across three depth layers and two size fractions [32], we analyzed the effects of phylogeny and niches on microbial co-occurrence networks across the levels of nodes, motifs, modules and the whole network, and further developed a novel mathematical model for simulating the self-organization of microbial networks based on a scale-free network [33,34]. We focused on three key questions: (1) How are direct and indirect interactions in microbial co-occurrence network constrained by both evolutionary and ecological processes? (2) What is the relative importance of evolutionary and ecological processes in determining co-occurrence network assemblages? (3) How do the characteristics of these two processes underpin co-occurrence network stability under environmental disturbance?

2. Materials and Methods

2.1. Datasets

We obtained meta data from the global ocean data of TARA project [32], and the samples were grouped into six subsets according to the three depth layers (SRF: surface water layer; DCM: deep chlorophyll maximum layer; MES, mesopelagic zone) and two size fractions (free-living: 0.22–3.0 μm; particle-associated: 0.80–5.0 μm) [32,35]. For the consistency of downstream analyses, the read depths of each sample of metagenomic 16S rRNA tags were rarefied to 10,000. Phylogenetic tree of representative OTUs was extracted from SILVA 128 database (QIIME release) [36] based on 97% similarity.

2.2. Network Construction, Motif and Latent Space Analyses

In each sub-dataset, OTUs observed in over a half of the samples were selected for network constructions. The correlations between OTUs were inferenced via SparCC tool [37,38,39] and carried out by ‘sparcc’ function (iter = 100, inner_iter = 100, th = 0.3) [38]. The correlations between OTUs that were no less than 0.65 were kept for network constructions. The memberships of nodes were identified by the fast-greedy modularity optimization method [23] via ‘cluster_fast_greedy’ function. The frequencies of each node for motif positions in undirected networks were counted using Simmons’s method, which was conducted via the ‘mcount’ function (six_node = TRUE) [40]. For each node in the network with the largest connectivity, the positions within the latent space in two dimensionalities and the probabilities of clustering memberships were calculated by the Variational Bayes Latent Position Cluster Model (VBLPC) [41]. To determine the optimal number of groups, we searched the minimum Bayes information criterion (BIC) with desired group scopes (Figure S1A). The latent space distances between OTUs were calculated as Euclidean distances. These network analyses were conducted in R version 3.6.1 with the packages of ‘SpiecEasi’ V1.0.7, ‘iGraph’ V1.2.4.2, ‘bmotif’ V1.0.0 and ‘VBLPCM’ V2.4.5 [38,40,41,42].

2.3. Statistical Analyses

Phylogenetic tree was transferred into pairwise distance matrix via ‘cophenetic’ function. The optimal habitat value, as a niche for each out, was calculated based on environmental variables [43]. Between-OTU niche differences were calculated as Euclidean distances, and between-OTU differences in topology properties of network motifs and module levels were calculated as Bray distance based on motif position frequencies and the probabilities of clustering memberships. We used multiple logistic regression in the distance matrices with 999 permutations [44] to quantify the effects of phylogeny and niches on network adjacent matrices (binary data). Similarly, we also used multiple linear regression in distance matrices to examine the regression coefficients of pairwise phylogeny and habitat distance against pairwise motif or module differences. These analyses were performed with packages of ‘base’, ‘vegan’ V2.5–6 [45] and ‘ecodist’ V 2.0.1 [46].

To estimate species preference in phylogeny and niches within networks, we proposed the novelty network-based metrics (Figure 1), that is, the mean neighbor phylogeny distance () and mean neighbor niche distance (), respectively. These two metrics were calculated by the following formula: and: where and denote phylogenetic and niche distance between species and its neighbor , and is the number of neighbors belonging to species . Further, we constructed 1000 random networks corresponding to empirical networks according to Erdos-Renyi model and generated a null distribution of and values. The standardized effect scores of and were given by: and:

Figure 1

The flowchart of the summarizing procedure for estimating mean neighbor phylogeny and niche distance (MNPD and MNND) at node level. First, the observed MNPD and MNND (obs.MNPD and obs.MNND) were quantified based on empirical networks. Second, the links of empirical networks were randomly shifting with the Erdos–Renyi model, and null distributions of MNPD and MNND (null.MNPD and null.MNND) were generated. Third, the standardized effect sizes of MNPD and MNND (ses.MNPD and ses.MNND) were calculated based on observed and null values. The absolute values of ses.MNPD (ses.MNND) larger than two show statistical significance (two sides, p < 0.05).

To further explore the design principles of ecological networks at the motif level, we calculated the preferences in phylogeny and niche of motif position () by: and: where is the frequency of the node in motif position (). As for the module level, the preferences in the phylogeny and niches of module were given by: and: where is the probability of node clustering into module .

We tested the linear relationships between the phylogenetic or niches distance and network path distance with ‘lm’ functions. Further, a random forest algorithm [47] was used to classify the network path distance based on species phylogeny and niche distance. Finally, we calculated the Pearson’s correlations between phylogeny or niches distance and latent space distance across network path distance. These analyses were performed with packages of ‘base’ and ‘randomForest’ V4.6–14.

2.4. Mathematical Mechanism Model

2.4.1. The dynamic of a Scale-Free Network in Ecology Community

Most of ecological networks were generally self-organized and exhibited scale-free property, such as food webs [48,49], species co-occurrence [20,50] and pollination networks [51,52]. In scale-free networks, the probability of a vertex with certain links (or degree) followed a power law distribution [33,34,53]: where is an attraction factor. A greater means the weaker ability of a node to attract other nodes [54]. The value of < 3.47 holds for almost all of empirical networks, and the most of these are range from 2 to 3 [53,55].

To explore the assemblage dynamic of microbial co-occurrence networks, we took evolutionary and ecological processes into account and extended the above formula (Figure 2A). The connectivity strength between species and is given by: where and are the weight of preferences for the similarity of phylogeny and niches to , respectively, and is the average degree of and . The parameters of and are the mutual attraction factors between and driven by phylogeny and niche similarity, as described by following formula: and: where is a constant representing the original attraction between species, and are the weights of the preference of phylogeny and niches similarity among species across network path distance . If , and reflect the direct phylogeny () and niche () attraction, respectively, originating from the preferences for phylogeny and niche similarity between and . If , and reflect the weight of the indirect phylogeny () and niche () attraction, respectively; attraction received from neighboring (); and the non-neighboring () vertexes with a distance of to or .

Figure 2

The simulation of self-organized assemblages in microbial networks. (A) The self-organized assemblages of microbial networks constrained by species phylogeny and niches. The connectivity c between species i and j is constrained by phylogeny (P) and niches (N). These two constraints could be directly affected by this pair of species (l = 0) and their neighbor (l = 1) or non-neighbor species (l > 1). (B) Characteristics of species pool and network initialization. In species pool, species i has preferences in phylogeny (Pp) or niche similarity (Np). The phylogeny and niche distance between species i and j are denoted by PD and ND, respectively. For network initialization, three pairs of species are randomly selected and linked. The initialized network is denoted by symmetrical adjacent matrix M (t = 0). (C) The simulation of self-organized assemblages in microbial networks. Symmetrical adjacent matrix M (t) was repeatedly updated according to network assembly procedures until termination events were triggered (See details in Supplementary Material) and the complete assembly of a pre-steady network. When environment heterogeneity changes (ΔEh), M (t) becomes involved in iterations until iteration termination and achieves a post-steady state. (D) The link dynamic in network assembly process. The red line is the number of network links across assembly time. The green and blue lines are the number of loss and gain links in assembly processes.

We assumed that and accorded with an exponential attenuation model: and: where and reflects the relative importance of phylogeny and niches, respectively. and is the attenuation rate of weight across network path distance , is the degree of environment heterogeneity. Then, is given by: and:

The mutual attractions between pair nodes driven by phylogeny () and niches () in network path distance are given by following formula: and: where and represent the direct mutual attractions of phylogeny and niches between and , respectively, given by: and: where and denote the preference of species in phylogeny and niche similarity, respectively. The phylogenetic and niche distances between and are represented by and , respectively. Considering the indirect effects on the co-occurrence of and , the phylogenetic and niche attractions for or from non-neighbor species of or with a distance of are given by Formulas (19) and (20), respectively, as follows: and: where . and ( is the aggregations of neighbor or not-neighbor nodes belonging to and with distance , respectively. The size of and ( are and , respectively.

2.4.2. The Specie Pool Construction

In order to obtain the aggregation of species for ecological network simulations, we generated species pools with regard to phylogeny and niches (Figure 2B). The community phylogeny was constructed by a random phylogenetic tree according to the Paradis’s algorithms [56], and the branch length was generated from uniform distribution . The niche distance between species and is given by: and: where is phylogeny signal in the optimal niche, which is a threshold of the phylogenetic distance for the linear relationship between phylogenetic and niche distance. The parameters of , and , are constant, and is noise sampled from the normal distribution . Finally, pairwise phylogenetic distance and niche distance are standardized by -score transformation.

2.4.3. Modeling of Network Assembly

At time , the species pool with species and a symmetrical adjacency matrix () with links is generated. Additionally, we then repeated the following five steps until there were no loss and gain links among nodes or the number of iterations reached the maximum value. The brief procedures are shown in Figure 2.

Step 1: At time , existing species in networks are likely to interact with new introduced species or other existing species that have no interaction with species . Regarding this, through Formula (25), we first calculated the pre-connectivity strength between connected node and its non-neighbor node with a pre-average degree of and : :

Step 2: To restrict the number of novel interactions generated at time , we define a threshold that is a quantile of , with a given probability and a minimal threshold . If > , node would prepare to connect node with probability through a binomial experiment. After that, a pre-connectivity adjacency matrix is created.

Step 3: Then, we generate a transitive adjacency matrix through:

Step 4: After the introduction of novel species or interactions, original network patterns are changed and some existing interactions between species are weakened and may disconnect. Thus, we calculated the connectivity strength matrix based on through Formula (25), and then defined a threshold of connectivity strength for disconnection. If < , the link between and would disconnect with a probability through the binomial experiment. After that, a pre-disconnected matrix was created.

Step5: Finally, we obtained the new adjacency matrix at time through:

2.4.4. Sensitive Analysis

We conducted a sensitivity analysis for the network assemblage model by Gaussian generalized linear model (GGLM) and random forest classification. First, the importance of phylogeny was set as twice as niche in the network assemblage, and the other of parameters varied within the given ranges (Table S1). Then, we varied the rest of parameters and generated 100 sets. For each parameter set, we sampled 100 replicates from modeling simulations. Finally, we defined the sensitivity of modeling outcomes to parameters as the regression coefficients of GGLM, where the dependent variables were the contributions of phylogeny and niche distance to network path distance given by random forest classification, and the independent variables were variable parameters.

To simulate the response of networks in the changes of environment heterogeneity, we examined the effects of the indirect interactions on the response of network dynamics in environment heterogeneity. In this model, environment heterogeneity is represented by , which alters the weight of the contributions of phylogeny and niches to the probability of species co-occurrence. After completion of pre-steady network with presetting (default ), was increased to 2 and adjacent matrix of was repeatedly updated. When the iterations achieved termination, we obtained a post-steady network . Finally, we calculated link turnover, including the percentages of link gain () and loss () by following formula: and: where represents the number of links in pre-steady network.

3. Results

Species phylogeny and niche similarity significantly correlated (p < 0.05) with, but showed differentiated contributions to, the marine microbial co-occurrence network structure at the motif, module and whole-network level (Figure 3A). Particularly, phylogeny made more contributions to network construction than niches at the motif level (Figure 3A; paired t-test, p = 0.07), while it was inconsistent at the module and whole-network level (Figure 3A; both p > 0.1). Specifically, at mesopelagic zone (MES), phylogeny, rather than niches, was a primary driver of particle-associated communities across three network levels (Figure 3A; p = 0.016), but not for free-living communities (Figure 3A). When the network path distance between species was considered, both phylogeny and niche distance positively correlated with network path distance up to around three, and then showed incongruent trends (Figure 3B). Such a phenomenon was further supported by the fact that phylogeny and niche distances were significantly (p < 0.05) positively related to network path distance within the same modules (Figure 3C), and these relationships were stronger than those across modules (Figure 3C). These results implied phylogeny and niche processes had strong constraints on network assemblage within a short range of network path distance, which may be particularly important for network module clustering.

Figure 3

The effects of species phylogeny and niches on marine microbial network constructions. (A) The relative contributions of species phylogeny (Phylo) and niche distance to network construction across motif, module and whole-network (whole) levels in free-living (Free) and particle-associated (Part) communities at surface water layer (SRF), deep chlorophyll maximum layer (DCM) and mesopelagic zone (MES), where pattern size and color intensity are the strength of relative contribution, and colors represent the positively (blue) or negatively (red) relative contributions. (B) The relationships between network path distance and phylogeny or niche distance at the whole network level, where color represent depth layers (red: SRF; green: DCM blue: MES), and line type represents community fractions (solid: Free; dashed: Part). The curves were fitted by generalized additive models. The shaded region indicates standard errors of the curves. (C) The linear relationships between network path distance and phylogeny or niche distance across modules (blue line) or within the same modules (red line), where line type is the statistical significance of linear regression analysis (solid: p < 0.05).

To determine whether species prefer to interact with phylogenetically close species or ecologically similar species, we developed novel indicators of mean neighbor phylogeny distance (MNPD)- and mean neighbor niche distance (MNND)-based ecological networks. Compared with the module level, nodes and motif positions showed the highest and lowest divergence, respectively, in ses.MNPD-ses.MNND plots (Figure 4A). Across six groups of microbial communities, over 95% of nodes showed a preference in phylogeny or niche-similar neighbors (ses.MNPD < 0 or ses.MNND < 0) and more than 65% of them showed a statistical significance (ses.MNPD < −2 or ses.MNND < −2; Figure 4B). For motif positions, all ses.MNPD and ses.MNND values were less than 0. This was especially true for the mesopelagic zone as all motif positions showed significant preferences in phylogeny or niche similarity: ses.MNPD or ses.MNND values < −2 (Figure 4B). These results indicate that phylogeny and niche processes have stronger constraints on network motifs or modules than on nodes. In addition, ses.MNPD had more significantly (p < 0.05) positive influences on node connectivity in particle-associated communities than in free-living ones, which were supported by the generalized linear model (Table 1).

Figure 4

The mean neighbor phylogeny and niche distance of species in microbial networks. (A) The ses.MNPD–ses.MNND plot for marine microbial communities across node, motif and module levels, where ses.MNPD or ses.MNND of motifs and modules were calculated according to the proportions of their node members. (B) The block proportions in ses.MNPD–ses.MNND plots across nodes, motifs and modules, where a ses.MNPD (or ses.MNND) larger than 2 indicates that one object has a significant preference in a more distantly related phylogeny (or niche) at same network level, while a ses.MNPD (or ses.MNND) less than −2 indicates one object has significant preference in a more closely related phylogeny (or niche).

Table 1

The influences of ses.MNPD and ses.MNND on node connectivity were analyzed by the generalized linear model. SRF: surface water layer; DCM: deep chlorophyll maximum layer; MES: mesopelagic zone.

Size Fraction	Layer	ses.MNPD	ses.MNND
Free-living	SRF	0.133	−0.0236
	DCM	0.2302	0.8291 *
	MES	−0.1677	−0.2277
Particle-associated	SRF	0.7122 ***	−0.0927
	DCM	0.3903 **	−0.0811
	MES	0.1315	−0.7282 **

***: p < 0.001; **: p < 0.01; *: p < 0.05.

Moreover, a random forest classification analysis [47] showed that phylogeny and niche distance affected both direct and indirect species interactions, that is, network path distances = 1 and >1, respectively (Figure 5). Regarding a large network path distance, the contributions of phylogeny and niches for the network path distance between two species generally increased and then decreased with a network path distance of around six. Such patterns showed a peak in the network path distances from two to five (Figure 5A). Notably, the contribution of phylogeny was around two times higher than that of niches in particle-associated communities at MES, and showed stronger long-lasting effects on indirect species interactions (Figure 5A). The highest importance of phylogeny and niches to network path distance classification were 0.286 and 0.247, respectively (Figure 5A).

Figure 5

The contribution of species phylogeny and niche distance to their direct and indirect interactions. (A) The relative contribution of species phylogeny and niche distance to classifications of network path distance estimated by random forest classification analysis. The colors of points or lines vary with phylogeny (orange) and niche (green). The line is a Loess fit based on the points. (B) The Pearson’s correlations between species phylogeny or niche distance and network latent space distance ignore the potential for modularity across the network path distance. The colors of bars or lines vary with phylogeny (orange) and niche (green). The lines are a Loess fit based on the values of bars.

To further assess the influence of phylogeny and niches on indirect species interactions, we calculated the latent space distance between species with a variational Bayes method [41] and estimated the Pearson’s coefficients between the distance matrices of latent space and phylogeny or niches across the network path distance. We found that species phylogeny or niche distance significantly correlated with the latent space distance (p < 0.05; Figure 5B), and these correlations generally decreased with increasing network path distance when ignoring modularity (Figure 5B). Interestingly, when including modularity, the constraints of phylogeny on species latent distance in ecological networks weakened. This is especially true for the particle-associated communities in the deep chlorophyll maximum layer (DCM; Figure S1B) and there were no significant (p > 0.05) correlations between species phylogenetic and latent space distances.

To disentangle the ecological significance of the observed indirect species interactions from the perspective of network stability, we established a mathematical model based on scale-free network and simulated the processes of microbial co-occurrence network self-organized assemblages (Figure 2). In a sensitive analysis of model parameters, our results showed that the distance of long-lasting effects l had a greater influence than the other parameters, such as community size S and attraction constant a, regarding the contribution of phylogeny and niche to indirect interactions, particularly in long-distance indirect interactions (Tables S2 and S3). In addition, the strength of the phylogeny signal in the niche optima of species significantly (p < 0.05) enhanced the contribution of niches to direct interactions in the community assemblages, primarily driven by phylogeny (Table S3).

We further simulated the response of the ecological network in environment heterogeneity (defined as continuous variable ), which could decease the weight of phylogeny preference on the probability of species co-occurrence [57]. After completeness of pre-steady network assemblage, was shifted from 0 to 20 and the post-steady network and was completed until the iterations of the matrix of connection strength terminated again (Figure 2C). To estimate network stability, we then calculated the link turnover [58], that is, the percentage of the gain and loss links. When species indirect interactions were considered in network construction, the ecological network had a better stability in more heterogeneous environments. Specifically, towards the long distance of long-lasting effects l, the percentage of gain (Figure 6A) and loss (Figure 6B) links reduced. Further, the influences of long-lasting effects on network stability in reducing link turnover were stronger in large, as opposed to small, communities (ANOVA p < 0.05). Moreover, the percentage of gain links showed unimodal patterns along the gradient of phylogenetic signals in optimal niches, with a peak at around 0.08 (Figure 6C). For the percentage of loss links, however, there was a minor variation until the phylogenetic signal increased up to around 0.01, and then it decreased rapidly (Figure 6D).

Figure 6

The response of network links in the changes in environmental heterogeneity. The gain links (A) and loss links (B) of the network across the distance of long-lasting effects (l) and community size (N; labeled with different colors) when environmental heterogeneity changes. The lines are a Loess fit based on 600 iteration points (6 × 100 replicates). The shaded regions indicate standard errors of the Loess curves. The gain links (C) and loss links (D) of the network across the phylogenetic signal in optimal niches (here, l = 2 and S = 50). The lines are fitted by generalized additive models (GAM) based on 10,000 iteration points (100 × 100 replicates). The shaded regions indicate standard errors of the curves of GAM.

4. Discussion

The evolutionary and ecological processes on direct [5,6,8,59,60] and indirect species interactions [1,61,62,63,64] attracted extensive attention; however, the relative importance of these two processes on ecological networks and its assemblage dynamics are less understood. Here, we clearly showed that species interactions were evolutionarily and ecologically constrained in marine microbial co-occurrence networks across the four levels of nodes, motifs, modules and whole networks. We further found that evolutionary processes revealed stronger long-lasting effects on indirect species interactions than ecological processes and dominated the network assemblage in spatially homogeneous but microscopically heterogeneous environments, such as particle-associated communities in MES. Finally, we simulated the assemblage dynamics of microbial co-occurrence networks with mathematically modeling and revealed that, regarding large network path distances, long-lasting effects improved network stability, as indicated by the decreasing link turnover.

We found that evolutionary and ecological processes generally constrained direct and indirect species interactions in marine microbial co-occurrence networks. For direct interactions, the findings were consistent with previous reports on predator–prey networks and competition networks [11,60]. For instance, over 80% of the observed realized and unrealized predator–prey interactions, involving 20 ground beetles and 115 prey species, could be predicted with phylogeny and functional trait information [11]. We further found that node connectivity in ecological networks was also associated with phylogeny and niches, implying that there was a phylogenetical conservatism [65] or niche specialization [66] in species interactions.

For indirect interactions, the findings on the strong effects of the two processes were consistent with previous reports of empirical networks among trophic levels, such as plant–animal-mutualistic [67] or predator–prey networks [1]. For example, with regard to the evolutionary processes in mutualistic networks, it has stronger constraints on consumer levels than resource levels and maintains network stability via positive within-guild indirect interactions between phylogenetically related species [67,68]. As for ecological processes, species density and function traits mediate indirect species coexistence and coextinction in plant–insect communities [1,62,69]. However, these reports show that indirect species interactions merely happen among three or four trophic levels, while we highlight that the two species indirectly connected within an ecological network could be generally constrained by phylogeny and niches, no matter whether they belong to different trophic levels or not.

Species interactions are fundamental for the assembly of ecological networks, and we expect that the assembly of the co-occurrence network is also driven by evolutionary and ecological processes. However, the ways that co-occurrence network assembly is driven by these two processes is a top-level issue, and could help us to essentially understand the self-organization of the microbial co-occurrence network. With integrations of empirical results, we proposed three reasonable principles of co-occurrence network assembly in marine microbial communities.

First, the phylogeny and niche are the primary drivers used to build the basic network blocks, such as motifs. For species, our null model analyses suggested that over 74% of species preferred to co-occur with phylogeny- or niche-related species, which was strongly supported by previous studies [20,26,28,70]. For motifs, all motif positions displayed the accordant preferences in phylogeny and niche similarity, which revealed that evolutionary and ecological processes had strong influences on building network motifs. It should be noted that, here we did not further consider that the effects of phylogeny and niches may vary regarding the type and strength of species interactions [10,12,71,72,73,74].

Second, both drivers had long-lasting effects on indirect species interactions within co-occurrence networks, and the strength of such effects decayed toward large network path distances. The long-lasting effects observed here are that the phylogeny and niche distance of pair species partly contribute to their indirect interactions, despite the large network path distance between them. Such findings agree with previous studies on food webs [12,63,75]. Furthermore, we revealed that the strength of long-lasting effects generally weakened with the increasing network path distance. Interestingly, the long-lasting effects of phylogeny and niches mostly occurred within a network path distance of around six, which may explain the phenomenon that the long-lasting effects were stronger within the same modules than across different modules. Overall, these long-lasting effects are short-range and could be the key mechanisms to enhance network modularity.

Third, phylogeny was more important than niches in determining direct or indirect species interactions when spatial homogeneity was greater. As marine depth increased and spatial homogeneity increased, the relative contributions of phylogeny to network constructions increased, while niches decreased (Figure 3A). For example, in the mesopelagic zone and surface water layer, phylogeny dominated the network assembly of particle-associated communities. However, free-living communities in the mesopelagic zone have no such phenomenon, which reveals that the microscale stability provided by particle may further strengthen the effects of phylogeny on the co-occurrence network assembly.

Based on these three principles, we proposed a mathematical model based on a scale-free network and simulated the self-organization of the co-occurrence network, which was key to exploring the network properties preceding empirical evidence. For example, regarding large network path distance, the long-lasting effects of both processes reduced link turnover and gave the network a strong resistance to environmental disturbances, indicating that the effects of indirect interactions on network stability may be more important than direct interactions [76]. This is particularly true for trophic networks with a low connectivity and high modularity [77] and for mutualistic networks with a nested architecture [78,79]. Moreover, the strong phylogenetic signals in optimal niches stabilized the successions of the co-occurrence network as expected when the determinants of network assembly shifted from evolutionary to ecological processes [80] (Figure 5B). However, for weaker phylogenetic signals, the percentage of gain links increased with the increasing phylogenetically signal. This might be caused by the increasing level of species connectivity with a higher strength than the minimum threshold of connections. This priori knowledge highlights that network stability should be further explored from design principles, except for network topologies such as size [81], connections [82] and modularity [83], because the latter are post analytical results and the outcomes of the former. Finally, the contributions of phylogeny and niches to direct and indirect species interactions had a great variety despite being from the same set of parameters, and this variety was greater for large network path distances (Figure S2). On the one hand, this supported multi-stability as the mechanism gave rise to different communities under the same environmental conditions [84]. On the other hand, it also indicated that stochastic processes may greatly affect the contributions of phylogeny and niches to network assembly in remotely indirect interactions.

Here, although we show a reasonable picture on how species phylogeny and niches shape microbial co-occurrence networks, three major caveats could be acknowledged regarding the analytical framework. First, microbial co-occurrence networks were inferenced by the correlation detection technique, which would lead to reasonable false positive rates [39], although we did not consider extremely rare species in network constructions. Meanwhile, specific ecological relationships, such as predation, competition, symbiosis and mutualism, could not be further accessed, and these ecological networks may have different and novel findings. Second, there are many correlation detection methods used to construct microbial networks, such as CoNet [85], RMT [86], MIC [87], LSA [88] and FlashWeave [89], which consistently produce different numbers and types of significant species correlation for the same data. Here, we used the SparCC tool for co-occurrence network constructions because it performed well in compositional data and data sparsity [39] and had significant edges with high correlations for analyses, which could greatly avoid noise and better than other tools. Third, some direct species interactions could be driven by indirect species interactions through other intermediary species, which could not be completely detected by the current methods and thus would decrease the robustness of results.

5. Conclusions

We found that direct and indirect species interactions in marine microbial co-occurrence networks were evolutionarily and ecologically constrained across four network levels of node, motif, module and the whole network, resulting in the small-world phenomenon of ecological networks. Compared to ecological processes, evolutionary processes dominated network assemblages in spatially homogeneous environments and showed stronger long-lasting effects on indirect species interactions. The contributions of two processes decreased toward long network path distance, but such long-lasting effects played a key role in stabilizing network link turnovers under environmental disturbances. These results provide a novel insight into the self-organization and modularization of marine microbial co-occurrence networks and highlight the importance of species phylogeny and niches in ecological network assemblages.