The eco-evolutionary assembly of complex communities with multiple interaction types

Domain: chemistry
Published: 2026-03-05
Editorial outcome: minor_revision

Abstract
Identifying the mechanisms that generate structure in complex ecological communities is fundamental for understanding their assembly, yet a comprehensive picture of how ecology and evolution combine to generate these patterns remains limited. We use an eco-evolutionary model of community assembly that incorporates interaction-driven population dynamics and evolutionary processes, including speciation and inheritance of interactions, to unveil the mechanisms generating and maintaining biodiversity in complex species interaction networks. Importantly, our model unpicks the effects of selection of interaction types from those of inheritance by comparing evolutionary assembly with invasion-based assembly under different combinations of interaction types. We find that a cost-benefit balance in accumulating interactions separates communities into two distinct types. Weakly beneficial interactions produce sparse, competition-dominated networks, whereas strongly beneficial interactions generate highly mutualistic, more connected communities. Mutualism, driven by both selection and inheritance, facilitates the emergence of large communities with increased complexity. Comparing model results with empirical patterns from microbial communities, we identify potential drivers of ecosystem assembly and characteristic interaction structures. Our results provide a classification system of complex ecosystems based on their composition of ecological interactions, thus generating testable hypotheses on the conditions under which different community types (mutualistic vs. competitive) might emerge.

Introduction
Unveiling and understanding the processes that generate and maintain biodiversity is a central focus of ecological research. Both evolutionary (e.g. speciation) and ecological (e.g. invasions) assembly mechanisms, and how they are shaped by biotic and abiotic factors such as species interactions and environmental conditions, have been shown to promote the complexity observed in natural communities [1–6]. Clarifying the relative importance of different mechanisms that foster biodiversity, while maintaining the stability of complex communities thus assembled, is fundamental to understand ecosystems [7, 8] and to inform future conservation strategies [9, 10].

Five decades ago, Robert May [11] demonstrated that ecological systems with a large number of species and high connectivity are less likely to be stable than simpler ones. This seminal finding sparked a decades-long search for the organisational features that enable natural communities to increase in complexity without compromising stability or persistence. May’s result referred to randomly connected communities (i.e., not subjected to a deterministic structure or assembly processes) represented by randomly sampled community matrices. Subsequent work has shown that particular network structures can alter the complexity-stability relationship [12], and that in empirical systems a large number of species often translates into low connectivity between them [3]. Furthermore, beyond structural patterns, the relative composition of interaction types, such as mutualism, competition, and consumer-resource interactions, has emerged as a key factor influencing the outcomes of community assembly and the stability of the resulting communities [13–17]. Only recently has theoretical research begun to explore how the interplay between different interaction types shapes the temporal assembly of complex communities, with recent studies revealing the key role of mutualistic interactions in fostering diversity (e.g. [18]). The next challenge is to uncover the mechanisms driving the emergence of different compositions of interaction types, the structural patterns of the resulting species interaction networks, and the influence of these patterns and processes on the diversity-connectivity relationships that underpin the complexity of ecological communities.

Evolutionarily, speciation is a major driver of community assembly. Speciation, acting on an existing pool of species and traits, together with subsequent selection, has been shown to recover patterns of biodiversity and network structure observed in complex model ecosystems [4, 19]. These findings suggest that speciation-driven assembly offers a key pathway for the emergence of complexity in ecological communities. Comparing speciation-driven evolutionary assembly with ecological assembly via repeated invasions helps clarify the distinct contributions of evolutionary processes to community structure [20]. Considering both speciation and invasion processes together, and their interplay, enables the investigation of two key mechanisms shaping community structure: selection of interaction types and inheritance of interactions. To date, theoretical studies of evolutionary community assembly have primarily focused on single-type species interaction networks such as food webs or mutualistic networks [4, 5, 19–25]. The combined effects of multiple interaction types have so far only been explored in a preliminary manner in evolutionary models of community assembly, with an emphasis on the study of stability [26], leaving their broader role in shaping ecological complexity largely unresolved.

Network structure not only determines biodiversity, complexity, and stability but can also characterise distinct types of ecological communities. These structural properties often emerge from different ecological and evolutionary assembly mechanisms, as observed across diverse ecosystems [27]. For instance, in microbial communities associated to phytoplankton cells in the phycosphere, community composition is shaped by host identity and metabolic profile [28, 29]. More generally, across host-associated microbiomes, microbial succession is thought to be driven by cross-feeding interactions between microbes and with the host [30]. The discovery of universal macroecological patterns of abundance and variability in microbial communities has motivated efforts to identify general network structures that can give rise to such patterns [31, 32]. Non-interactive population dynamic models have successfully reproduced many of these patterns and have been proposed as null models for microbial community structure [31, 33, 34]. However, such models fail to capture non-random correlations in species abundances observed in these communities, a caveat that has been addressed by adding sparse and weak interactions between community members [35, 36]. Together, these findings highlight the critical role of the structure of ecological interactions in shaping microbial communities and the need to identify the assembly mechanisms that give rise to such structural organisation. Similar challenges are observed across ecosystems types and scales. For example, in mutualistic networks comprising interactions between plants and their animal pollinators, characteristic patterns of network modularity and nestedness have been suggested to arise through adaptation and niche optimisation [25] or speciation [4]. A holistic eco-evolutionary understanding of how these processes and their interplay can shed further light on ecological structure across scales.

Here, we develop a deeper mechanistic understanding of community assembly by investigating the evolutionary emergence of complex ecological networks involving multiple interaction types. We use a dynamical model that we previously developed [26] in which mutualistic, consumer-resource, and competitive interactions shape the outcomes of speciation and invasion events driving community assembly and complexity. Interactions between species determine how they influence one another, potentially triggering extinction cascades. Over time, successive assembly events and ecological filtering give rise to stable interaction network structures. We find that the strength of ecological interaction benefits defines a threshold that separates two distinct community types: one dominated by competition, and the other by mutualism. Moreover, mutualism and speciation jointly promote increasing network complexity and species diversification. Nonetheless, species richness can increase independently of mutualism and speciation if enough invasion events occur and network connecitivity is low. By analysing longitudinal data from microbial communities [31, 35], we link macroecological patterns to underlying assembly mechanisms. We hypothesise that strong benefits from ecological interactions can favour the emergence of highly interactive microbial communities that display patterns of abundance correlations characteristic of these systems. By linking eco-evolutionary assembly mechanisms to emergent community structure and observed macroecological features, our results advance understanding of the processes that shape complexity in natural communities.

Results
We analyse a step-wise community assembly model in which new species are introduced sequentially to the existing community via two distinct mechanisms: (1) evolutionary speciation and (2) invasion. We start the assembly process with small communities comprising five non-interacting species.

Assembly by evolutionary speciation: In this scenario, new species arise as a variation of a randomly selected parent species and inherits its ecological interactions with limited fidelity. Inheritance with changes in species interactions represents mutations in foraging, morphological, or physiological traits that biologically determine interactions with other species (e.g. feeding modes, body morphology, metabolic pathways, signaling mechanisms, etc.). A parameter \Delta defines the maximum number of interactions by which the offspring may differ from the parent species. Smaller values of \Delta correspond to more constrained inheritance (i.e. greater similarity to the parent), a regime we term strong inheritance. For each speciation event, all interactions of parent species are copied into the offspring, although not faithfully: an integer d, representing the total number of interaction changes, is then drawn uniformly from the interval [1,\Delta]. Interactions are randomly removed and/or added from the offspring until the total number of modifications, relative to the parent, equals d. For a given d, the number of interactions created is uniformly chosen from [1,d] and the number of interactions destroyed is d minus that number. If there are not enough interactions inherited (or no interactions, as in the initial state), then new interactions are created whenever an inexistent one would be destroyed (and vice-versa). Speciation is thus a flexible process, with the average degree of inheritance controlled by \Delta: larger values lead to weaker average inheritance. If a parent species has no interactions (i.e. an isolated species), at least one interaction must be added (d \ge 1) to ensure variability.

Assembly by invasion: We then define a purely ecological scenario in which new species are introduced without any inherited interactions. Instead, each invader is assigned a connectivity value drawn uniformly from the interval [\rho_1, \rho_2]. This value defines the probability of forming an interaction with each existing species in the community. As invaders differ in their connectivity, community-level connectivity can shift over time due to ecological filtering / selection. We set \rho_1 = 0.01 and \rho_2 = 0.5 to span a range comparable to connectivity values observed in speciation scenarios.

New species are introduced once community dynamics reach equilibrium, i.e., when the dynamics governed by Eq. (1) have stabilised. Species addition can destabilise the system, driving existing species below an extinction threshold and leading to their removal from the community. The addition of a new thriving species (i.e. that does not go extinct) constitutes an assembly event (Fig 1a). After a new species is introduced in an assembly event, we assume no rewiring of its interactions, only the possible addition of interactions with new species included in subsequent events. Species abundances follow ecological dynamics governed by the system of differential equations in Eq. (1), where x_i is the abundance of species i in a community of S species, governed by several mechanistic drivers (Fig 1b). The non-negative coefficients p_{ij}, m_{ij}, and c_{ij} represent the strengths of consumer-resource (+/-), mutualistic (+/+), and competitive (-/-) interactions, respectively, and are nonzero only when there is an interaction of the corresponding type between i and j. These interactions are interspecific, thus p_{ii} = m_{ii} = c_{ii} = 0. These coefficients determine the strength of the effect of species j on the density of species i. For consumer-resource interactions, we have p^+_{ij} if i is a consumer and p^-_{ij} if i is a resource.

(1)

The variable n_i + m_i measures the number of interactions in which species i is a mutualist or a consumer, used to consider the cost of every positive interaction. The sums over j and k run for all S species. We used Type II functional responses to represent intake saturation in both mutualistic and consumer interactions. Note the subtle difference in the denominators of resource and consumer interaction terms, summing respectively over p^+_{jk} and p^+_{ik}, because the sum gathers all the resources of a consumer (i.e. all sources). Each species i has a positive intrinsic growth rate r_i drawn from a half-normal distribution (\mu_r, \sigma_r^2). Intraspecific competition strength s_i is chosen such that the resulting carrying capacities r_i/s_i follow a lognormal distribution. New interactions introduced via speciation or invasion are assigned a strength drawn from a half-normal distribution (0, \sigma^2), where \sigma sets the level of interaction strength (the choice of interaction type depends on the scenario). In the case of consumer-resource interactions, we impose a directional constraint on interaction strengths to reflect energetic inefficiency. Specifically, we assume that the number of resource individuals removed through consumption always exceeds the number of new consumer individuals produced as a consequence of that consumption. Accordingly, when interaction strengths are randomly sampled, if p^+_{ij} > p^-_{ji}, we reset p^+_{ij} so that p^+_{ij} = p^-_{ji}. This constraint prevents unrealistically efficient trophic conversion.

To capture saturation of resource intake, we assign two parameters: h_p for consumers and h_m for mutualists, representing their average handling times. Every positive interaction (whether a consumer-resource link in which species i is the consumer, or a mutualistic link in which i benefits) incurs a harvesting cost \delta, representing the energetic and physiological investment required to acquire and maintain benefits from a given interaction. Such costs arise from the allocation of finite resources to interaction-related traits. For simplicity and interpretability, we assume a constant per-interaction cost, understood as a mean-field approximation capturing the average energetic burden of maintaining a positive interaction. This cost is subtracted from r_i for each beneficial interaction, reflecting a trade-off: as a species accumulates positive links, it becomes increasingly reliant on its partners and may ultimately depend on them for survival (i.e., obligate mutualism or increasing trophic dependence on external resources for consumers, with reduced viability in their absence). This choice reflects the allocation of energy in either interaction traits or survival and reproduction. The balance between the interaction strength parameter \sigma and the per-link cost \delta therefore determines the net benefit of positive interactions.

A separation between competitive and mutualistic communities

The interplay between the strength and the cost of interactions creates a threshold that separates competitive from mutualistic communities. After 500 evolutionary assembly steps, we observed a threshold in the benefit of ecological interactions that separates two distinct community types (Fig 2). When \sigma is low and \delta is high, resulting in weak interaction benefits, communities are dominated by competitive interactions (Fig 2c). We refer to this regime as Type 1 (T1). Conversely, when \sigma is high and \delta is low, producing strong benefits from ecological interactions, mutualisms predominate (Fig 2a). We refer to this regime as Type 2 (T2). Consumer-resource interactions reach their highest proportion at the interface between T1 and T2, corresponding to intermediate levels of interaction benefit (Fig 2b).

Species richness (i.e. the number of species in the network) decreases with increasing harvesting cost \delta (Fig 2d), while network connectance (C = L/S^2 on the binary directed network, a measure of connectivity) reflects the threshold between community types, being higher in T2 communities (Fig 2e). As a result, community complexity, defined as the product of species richness and connectance, substantially increases in mutualism-dominated T2 communities (Fig 2f). Additionally, we found that the intensity of intraspecific competition shifts the position of the threshold (Supplementary Fig S2), with higher values shifting the regime from T2 to T1, indicating that the benefits of interactions must be evaluated relative to the strength of self-regulation. The existence of this threshold highlights interaction benefits as key drivers of community assembly, shaping the composition of interaction types and mediating the emergence of complexity.

A simple mean-field argument can clarify the origin of the threshold. Consider a community structured by mutualistic and competitive interactions. The potential advantage of a species acquiring a new mutualistic link can be expressed by the inequality

(2)

where x is the average abundance of a species, \mu is the mean interaction strength, k_1 is the mean number of mutualistic links per species, and \delta is the cost per mutualistic interaction. This inequality formalises the intuitive idea that a new mutualistic interaction is selectively favourable only if the associated benefit exceeds its cost. We can further approximate the equilibrium abundance of a species using a mean-field description:

(3)

where k_2 represents the strength of interspecific competition and \gamma the mean intraspecific competition. Here, the focal species experiences the community as an average background, reducing the dynamics to a single uncoupled equation for x. By combining these two expressions, we can identify the parameter region in which the inequality is satisfied, given that x > 0 and x satisfies the equilibrium equation. This analysis reveals the existence of a threshold in k_1, which shifts toward lower values depending on the parameter set. In particular, higher costs \delta for a fixed small value of k_2 cause the inequality to fail at smaller k_1, highlighting a clear benefit-to-cost trade-off associated with acquiring additional mutualistic links (Supplementary Fig S1). Intraspecific competition, being central in determining the equilibrium density x, also impacts the threshold.

Selection and inheritance drive interaction type composition

Having established two distinct community types driven by interaction benefits, we next examine the role of inheritance in shaping community structure. To do so, it is necessary to distinguish inheritance effects from ecological filtering driven by the selection of interaction types that characterises the compositions of T1 and T2 communities. Both interaction-type selection and inheritance contribute to assembly under evolution, but only selection is present in assembly by invasion. By comparing evolutionary with invasion-driven assembly both considering the selection of interaction types, we can isolate the effects of inheritance.

To this end, we designed four assembly scenarios. These scenarios combine two community types (T1 and T2), two inheritance regimes (weak and strong, for evolution), and two levels of variability of interaction types (low and high, for invasion), yielding a total of eight scenarios (Fig 3a). The key distinction among these scenarios lies in how new species acquire their ecological interactions during assembly events.

Evo. (\Delta=5): evolution with strong inheritance. A new species arises via speciation from a randomly selected parent species. After copying the parent’s interactions, up to five additions or deletions are allowed (\Delta=5), resulting in high inheritance fidelity and close similarity between parent and offspring. The types of newly added interactions are determined by choosing the signs of each directed interaction (positive or negative) with equal probability.

Evo. (\Delta = 30): as above, but with \Delta = 30, offspring may differ by up to thirty interactions, producing low inheritance fidelity and greater divergence from the parent species. Intermediate evolutionary scenarios (\Delta = 14 and \Delta = 22) yielded a similar behaviour to this one (Supplementary Figs S3 and S4, Supplementary Note 1).

Inv. Rand.: invasion with fixed probabilities for each interaction type. Each new species is introduced with independently assigned interactions, where each link is randomly assigned to one of the four available interaction types (mutualist, competitor, consumer, or resource). Due to the random nature of this process, this will yield a probability of 0.25 for each of the interaction types. Because these probabilities are fixed, the proportion of interaction types shows low variability across invaders, limiting the scope for ecological filtering (i.e. selection cannot select species with larger proportions of specific interaction types because they do not exist). This scenario serves as a baseline and mirrors previous models of invasion-based assembly with mixed interactions [18].

Inv. Prop.: invasion with variable proportions of interaction types. For each new invader, the overall proportions of the four interaction types are randomly drawn from a uniform distribution. Interactions are then assigned accordingly. This introduces high variability in the species-level composition of interaction types, allowing ecological dynamics to selectively shape (via selection) the structure of community-level interactions composition.

Community types (T1 and T2): for every scenario, we chose particular points (\delta, \sigma) to represent T1 and T2 community regions, (\delta = 2.5 * 10^{-2}, \sigma = 0.05) for T1 (fifth column from the left and second row from the bottom in Fig 2) and (\delta = 1.0 * 10^{-2}, \sigma = 0.2) for T2 (second column from the left and fifth row from the bottom in Fig 2). We tested different points within the same regions and all yielded similar analyses (Supplementary Figs S5-S8, Supplementary Note 1).

After 500 assembly events and at ecological equilibrium, the eight scenarios produced a broad range of interaction-type compositions (Fig 3b). Across all scenarios, the threshold between T1 and T2 communities was preserved: competitive interactions dominated in T1, while mutualistic interactions predominated in T2. However, inheritance and interaction-type selection each independently amplified the divergence between community types. Strong inheritance (\Delta = 5) led to the clearest predominance of the corresponding selected interaction type (consumer-resource or mutualistic), as well as the largest variability in the composition of interaction types. These results underscore the importance of species-level inheritance in shaping community-level structure and in reinforcing the effects of selection across community types.

Evolution fosters diversity and complexity in assembled networks

To identify the drivers of community complexity, we recorded the trajectories of each assembly scenario across 500 assembly events, evaluated at ecological equilibrium (Fig 4a,c). The increase in complexity (S * C) was strongly dependent on community type, with T2 communities reaching substantially higher values across all scenarios. Across both community types however, both selection and inheritance independently contributed to increases in complexity. In T1 communities, Evo. (\Delta=30) and Inv. Rand. exhibited similar trajectories, due to the low connectance inhibiting selection from effectively influencing the outcome of assembly when \Delta is high. By contrast, Evo. (\Delta = 5) and Inv. Prop. communities reached higher average values of connectance and richness (Fig 4b), leading to a more sustained increase in complexity following a rapid initial increase (Fig 4a). Although species richness increased substantially in all communities, connectance remained consistently low. This suggests that when interaction benefits are low, resulting in competitive T1 communities, assembly tends to produce sparsely connected networks with limited complexity.

In T2 communities, the Inv. Rand. scenario produced communities with the lowest complexity, primarily due to reduced species richness (Fig 4c,d). In contrast, Inv. Prop. resulted in higher complexity due to increased species richness, despite exhibiting lower average connectance (Fig 4c,d). Evolutionary scenarios had a stronger impact on assembly outcomes than invasion. Both weak and strong inheritance produced similarly high levels of complexity, the highest among all scenarios (Fig 4b), but through distinct structural pathways. Evo. (\Delta = 5) produced communities with higher richness and lower connectance, while Evo. (\Delta = 30) yielded the opposite pattern (Fig 4d). Thus, the strength of inheritance acts as a mechanism shaping alternative routes to complexity. Moreover, evolutionary scenarios resulted in greater heterogeneity of communities across the richness vs. connectance space, compared to invasion (Fig 4d). Overall, high interaction benefits that promote mutualism are also consistently associated with increased community complexity [26]. Both the selection of interaction types and inheritance contributed to the predominance of mutualism and the resulting emergence of complexity.

Our results suggest that the effectiveness of inheritance in driving the emergence of complexity depends on community type, and that the degree of inheritance determines the specific pathway through which complexity arises. T2 communities exhibit clearer distinctions between assembly scenarios, highlighting the stronger influence of both selection and inheritance in these mutualism-dominated systems.

We extended our analysis to explore faster evolutionary assembly events, out of equilibrium (Supplementary Figs S9-S12, Supplementary Note 1), mixed assembly events in different proportions of evolution and invasion (Supplementary Figs S13 and S14), and varying levels of intrinsic growth rate (Supplementary Figs S15-S18, Supplementary Note 1) and intraspecific competition (Supplementary Figs S19-S22, Supplementary Note 1). All additional results contribute to corroborating the robustness of our conclusions. Moreover, these results show that evolution with high inheritance is more resilient than other scenarios in assembling a community that thrives under stricter conditions (Supplementary Figs S15, S17, and S21).

To further examine how interaction architecture (i.e. the structure of who interacts with whom) emerges during community assembly, independently of interaction sign, direction, or magnitude, we analysed the topological degree entropy (a measure of the evenness in the number of links across species) and network modularity of assembled communities. For this analysis, community networks were represented in an undirected and unweighted form, retaining only the presence or absence of interactions. This choice isolates structural properties of interaction architecture that are robust to dynamical fluctuations in interaction strengths, while necessarily discarding information on interaction type, directionality, and energetic asymmetries. We quantified the relative increase in structure using a score metric of relative increase over random networks, comparing each observed value Z to the average of 50 randomised communities Z_r with matching complexity. The score was defined as Z_s = (Z - Z_r)/Z_r (see Methods for details). In T1 communities, strong inheritance (\Delta = 5) produced a higher relative increase in both degree entropy and modularity. In T2 communities, both inheritance and selection contributed to increased entropy, but only inheritance led to a substantial increase in modularity (Supplementary Fig S23). Weighted metrics yield equivalent results for T2 communities, but T1 communities no longer exhibit an increase in degree entropy (Supplementary Fig S24). These findings suggest that inheritance promotes the emergence of more modular networks with heterogeneous distributions of the number of interactions across species, both organisational hallmarks of ecological networks [37].

Assembled communities recapitulate microbiome organisation

We compared our model results to empirical data from microbial communities in an attempt to identify the mechanisms, assembly processes, or community types that best approximate observed macroecological patterns in these ecosystems. We focused on patterns of relative abundance distributions found in both cross-sectional and time-series microbial community data using data from the human microbiome across different body sites (palm, gut, and mouth) and comprising 3,121 microbial types found in n=1,341 samples [31]. As shown elsewhere [31, 35], these data exhibit consistent statistical patterns in the standardised log-mean abundance distribution (MAD) and in pairwise abundance correlations (Fig 5). Importantly, deviations from random expectations in pairwise abundance correlations have been linked, in part, to ecological interactions [35]. To generate data from our model outputs comparable to those from microbial communities, we ran stochastic simulations of communities assembled from equilibrium after 500 events (see Methods for details). To quantify the comparison between correlation distributions, we used the Wasserstein distance, which measures the total displacement of probability mass required to completely match two distributions.

We used the Inv. Rand. scenario, in which interactions are assigned randomly, as a baseline for comparison. Due to sparsity of interactions, all T1 scenarios produced similar relative abundance and pairwise correlation patterns (Supplementary Fig S25, Supplementary Note 2), as illustrated in Fig 5a,b for the Inv. Rand. baseline. This suggests that, in T1 communities, variation in assembly mechanisms does not translate into differences in abundance patterns, which may instead be dominated by environmental stochasticity. Moreover, none of the T1 scenarios represented well the empirical MAD or pairwise correlation patterns.

In contrast to T1, T2 communities yielded a closer resemblance to the empirical patterns observed in microbial communities. The Inv. Rand. baseline scenario approximates the observable range of the microbial mean abundance distribution (MAD), though it fails to capture the upper tail of highly abundant taxa (Fig 5c). It also reproduces the overall shape of pairwise abundance correlations (Fig 5d). Weak evolutionary inheritance (Evo. \Delta = 30; Fig 5e,f) and the Inv. Prop. scenario (Supplementary Fig S26) yield an even closer match to empirical patterns. These scenarios better capture the higher density region of the MAD (Fig 5e) and approximate the empirical distribution of negative pairwise correlations, although they underestimate the frequency of positive correlations (Fig 5f).

Based on the Wasserstein distances (Supplementary Fig S27) for the correlation pattern, T2 communities provide a better match to the empirical data than T1 communities. In T2 scenarios, the simulated distributions are closer to the empirical ones, and also closer to those obtained from the corresponding simulations without interactions (i.e. logistic simulations of the same communities). In contrast, T1 communities produce simulated distributions that deviate more strongly from the empirical data than their non-interacting counterparts, indicating that interactions in T1 communities systematically worsen the agreement with observations. Thus, whereas interactions in T2 communities improve empirical concordance, interactions in T1 communities drive simulations further away from empirical patterns. Extending this analysis to additional scenarios explored in the Supplementary Material yielded consistent results: good agreement was recovered for T2 simulations with mixed evolution–invasion assembly (Supplementary Fig S28) and under reduced intraspecific competition (Supplementary Fig S29), while no T1 scenario produced a comparably good match. Therefore, our analysis identifies T2 communities as better representations of microbial communities, highlighting strong interaction benefits as a potential mechanism underlying their macroecological abundance patterns.

Methods
Simulations

All simulations begin with five species with no initial interactions. Species abundances at initialization are uniformly drawn from the interval [0, 0.02]. A species is considered extinct and removed from the system if its abundance falls below a threshold of x_{ext} = 10^{-6} (we verified that results are qualitatively equivalent for x_{ext} = 10^{-5} and x_{ext} = 10^{-7} and thus not sensitive to a particular choice). Each simulation proceeds for 500 assembly events. Assembly occurs only after the system reaches ecological equilibrium, defined as a relative change in all species abundances of less than 0.01% over ten time steps of size 0.01. Intrinsic growth rates r_i are sampled from a normal distribution (\mu_r, \sigma_r^2) with \mu_r = \sigma_r = 0.1. Intraspecific competition coefficients s_i are determined such that carrying capacities K_i = r_i/s_i follow a lognormal distribution with underlying normal parameters \mu = 2.2, \sigma = 0.5; specifically, s_i \sim r_i / logN(2.2, 0.5^2). In evolutionary scenarios, r_i and s_i are not inherited (i.e. are sampled from their distributions, regardless of the parent), reflecting their dependence on species-specific environmental conditions and focusing inheritance on interaction structure only. Simulations use Type II functional response on positive interactions with constant handling times set as h_p = h_m = 0.1. When new interactions are introduced, their strengths are drawn from a half-normal distribution: (0, \sigma^2), where \sigma controls the typical strength of interactions. In evolution scenarios, inherited interaction strengths are subject to small stochastic variation, given by w' = w(1 + N(0, [0.05w]^2)), where w is the parent’s strength and w' is the offspring’s (interactions do not change type with inheritance). After the first 20 assembly events, all isolated species (i.e. those with no interactions) are removed after each subsequent event. Numerical integration of ecological dynamics (Eq1) is performed using the odeint function from the SciPy library in Python [51].

Analysis

For each combination of \delta = {0.005, 0.01, 0.015, 0.02, 0.025, 0.03} and \sigma = {0.01, 0.05, 0.1, 0.15, 0.2, 0.25} in the heatmaps, we ran 15 independent simulations using the strong inheritance evolution model (\Delta=5), each consisting of 500 assembly events. Final community properties were averaged across these replicates and used to generate the heatmap data. To analyse the temporal dynamics of community composition (Fig 4), we ran n=20 replicate simulations per scenario. Variables were recorded every 50 assembly events, and both the mean and standard error (defined as the standard deviation divided by the square root of the number of replicates) were computed and shown in the plots. Except for the heatmaps, all other results were based on these n=20 replicates per scenario. Type 1 and Type 2 communities were defined by parameter combinations (\delta = 0.025, \sigma = 0.05) and (\delta = 0.01, \sigma = 0.2), respectively (which are contained in the heatmaps).

Network metrics

We evaluated two structural properties of the assembled ecological networks: degree entropy and modularity. Degree entropy quantifies the heterogeneity of the network’s degree distribution P(k), where k denotes the number of interactions (degree) of a node (species). We computed the entropy using the unweighted network (i.e., considering only presence-absence of links), as H = - \sum P(k)ln(P(k)). Modularity is a standard network metric that quantifies the extent to which nodes form communities, i.e., groups of species that are more densely connected internally than to other parts of the network [52]. We calculated modularity from the unweighted interaction network using Mod = \sum_c [L_c/m - (k_c/2m)^2], where c indexes each community (or module), L_c is the number of intra-community links, k_c is the sum of degrees of nodes in community c, and m is the total number of edges in the network. The resolution parameter \epsilon was set to 1. All network analyses were performed using the Python package NetworkX [53], and community detection was carried out using the Louvain algorithm [54]. Effective increase: To facilitate comparison across scenarios, topological network metrics (degree entropy and modularity) were expressed as their effective increase relative to corresponding random networks. Specifically, for a given metric Z, we computed its effective increase Z_s as Z_s = (Z - Z_r)/Z_r, where Z is the observed value in the assembled network, and Z_r is the mean value obtained from 50 Erd\H{o}s-R\'{e}nyi random networks with the same number of species S and connectance C. This measure represents the percentage increase of the observed network property relative to what would be expected by chance. For the weighted comparison (Supplementary Fig S24), we maintained information about interaction strength and summed the absolute strength (m, p^+, p^-, c) of the pair of directed interactions between two species to compose the weighted link.

Data analysis

We utilised publicly available data from EBI Metagenomics: longitudinal human microbiome data (feces, left/right palm, tongue; subjects F4 and M3; Study MGYS00002184, ERP021896) from the EBI MGnify platform (project ERP021896) [55, 56]. This is the same longitudinal dataset used by Grilli (2020) [31] to investigate macroecological patterns. We retained only samples containing at least 10,000 sequencing reads and computed species’ relative abundances accordingly. The 6 filtered data groups utilised in the analysis contained n=131 (feces F4), n=134 (l-palm F4), n=334 (feces M3), n=143 (l-palm M3), n=268 (r-palm M3), and n=331 (tongue M3) samples. Two groups (tongue F4 and r-palm F4) were excluded due to their small resulting sample size. To generate comparable simulation data, we reused the same replicate simulations from Figs 3 and 4, starting from the final equilibrium state reached after 500 assembly events. We modelled species abundance fluctuations using stochastic differential equations with environmental noise. Specifically, we added independent Gaussian noise to each species, scaled proportionally to its abundance, with a standard deviation of 0.1. Integration was performed using the Euler-Maruyama method, and abundances were recorded at evenly spaced, non-overlapping time intervals. We computed Wasserstein distances using the function wasserstein_distance from scipy.stats in Python. The simulations without interactions (logistic) used for comparison were done by integrating stochastic perturbations with only the intrinsic growth rates and intraspecific competition in the equation for the assembled communities. We waited until a new equilibrium was reached before sampling the timeseries.

Data availability

No new experimental or field data were collected. Publicly available datasets from the EBI MGnify platform were used for model comparison, specifically longitudinal human microbiome data from subjects F4 and M3 in MGYS00002184 [https://www.ebi.ac.uk/metagenomics/projects/MGYS00002184]. Raw OTU data were filtered and transformed as described in the Methods, and the original tables and process are available in the Zenodo repository [https://doi.org/10.5281/zenodo.15621027].

Code availability

Files generated and analysed in this study and computer code developed to implement the model and execute the numerical simulations are available in the accompanying Zenodo repository [https://doi.org/10.5281/zenodo.15621027] and on GitHub: [https://github.com/computational-ecology-lab/eco-evo-network-assembly].

All simulations and data analyses were performed in Python (v3.11.11). The computational workflow relied on NumPy (v2.2.6) and SciPy (v1.15.2) for numerical computations, pandas (v2.2.3) for data handling, and Matplotlib (v3.9.1) and seaborn (v0.13.2) for visualisation. Additional specialised visualisations were generated using python-ternary (v1.0.8).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript is quite interesting and deals with a noteworthy topic in theoretical ecology. It shows two different types of communities emerging from assembly under different parameter regimes (interaction strength vs. harvesting cost), and then these two types are further investigated under different evolutionary inheritance (weak and strong) or variability of interaction types during invasions (low and high). This leads to 8 scenarios, for which complexity is examined in terms of species richness, connectance, and further network metrics. Relative abundance distributions are also compared to empirical observations on microbiomes.\nI have few suggestions for clarifications and improvements, listed below. Further details are also attached.\nHarvesting cost delta: this seems a relevant parameter in your model and investigation, but it is not very well explained or justified. Why is it needed? Is there any empirical or theoretical foundation for it? Why is it constant? The name itself is also not very informative.\nNew interactions for both evolutionary and invasion are extracted from a half-normal. The invasion scenarios Rand. and Prop. specify how the types of interactions are chosen, but it is unclear for the evolutionary scenarios.\nIntroduction of new species happens at equilibrium, are results robust if mutations and/or invasions were not rare?\nEvolution and invasion are treated separately, what results would emerge if both were acting at the same time?\nAmong the network metrics, nestedness is not mentioned, but often used and discussed in this type of research.\nAlso, unweighted metrics are used, are results robust to the inclusion of weights in network metrics?\nComparisons with empirical data are qualitative, but the similarity between simulations and data should rather be quantified formally.\nI hope these suggestions can help improve the manuscript.", "reviewer_2": "The eco‑evolutionary assembly of complex communities with multiple interaction types\nIn this study, the authors develop an eco‑evolutionary community assembly model that integrates three interaction types (competition, consumer–resource, and mutualism) within a single dynamical framework. Communities assemble via sequential introduction of species, either through speciation with inheritance of interaction structure or via invasion from a regional pool. A cost–benefit trade‑off is imposed on each positive interaction (mutualistic or consumer) via a per‑link harvesting cost and a saturating Type II functional response for the benefits. By exploring the plane of interaction strength (σ) and cost (δ), the authors identify a threshold separating two regimes: a competition‑dominated, sparsely connected “Type 1” community and a mutualism‑dominated, more connected “Type 2” community.\nThe paper then compares eight assembly scenarios that combine the two regimes (T1/T2), strong vs weak inheritance (Δ small vs large), and invasion with either fixed or variable interaction‑type composition. This is used to separate the effects of selection of interaction types from inheritance of interactions on emergent complexity (richness × connectance), degree entropy, and modularity. Finally, the authors compare time series generated from assembled communities to empirical macroecological patterns from human microbiome datasets, focusing on mean abundance distributions and pairwise abundance correlations.\nUpon review, I am convinced that combining multiple interaction types with explicit inheritance in a community‑assembly framework is timely and potentially impactful, especially given the connection to microbial macroecology. Generally, the model is well motivated and the results are interesting. However, I believe the manuscript can be substantially strengthened both conceptually and in clarity if the following major and minor points are addressed.\nMajor points\n1. The dynamical system is central to the paper, yet Eq. (1) is currently difficult to parse. As typeset, the different sums and terms (competition, consumer–resource, mutualism, costs) are condensed in a way that is hard to reconstruct, and some notation (e.g. which sums appear in the denominators, the meaning of “Σ_{p+,m} δ”, etc.) is not fully explained.\nI suggest:\n• Rewriting Eq. (1) in a more modular way (e.g. grouping explicitly into “self‑regulation”, “loss to consumers”, “gain from resources”, “gain from mutualists”, “harvesting cost” terms), possibly split over several lines.\n• Explicitly defining all summation indices: over which sets (consumers, resources, mutualists) each sum runs, and which interactions contribute to the denominators of the Type II functional responses.\n• Clarifying in the main text how consumer–resource asymmetry is handled (the rule p_{ij}⁺ ≤ p_{ji}⁻ and the enforced symmetry when this fails). This is an important modelling choice and deserves a brief justification, not only a note in Methods.\n• Stating explicitly whether the network is treated as directed and signed in the dynamics but undirected/unsigned when computing network metrics (see point 5).\n2. The authors show a sharp transition in interaction‑type composition, connectance, and “complexity” as σ and δ vary, and then select a single parameter pair in each region as “Type 1” and “Type 2”. Supplementary Fig. S1 further demonstrates that rescaling intraspecific competition shifts communities across this threshold, indicating that it reflects the balance of interaction benefits vs self‑regulation, rather than σ and δ alone.\nI suggest:\n• Emphasising that T1 and T2 are regimes or “community types” defined by this relative balance, not single parameter points. The selected (σ, δ) pairs should be clearly described as representative.\n• Providing at least a brief robustness check (in the SI) showing that key results for the assembly scenarios (e.g. the roles of inheritance and selection in Fig. 2–4) are qualitatively stable under modest changes in σ and δ within each regime.\n• Considering a simple heuristic or mean‑field argument for the location of the threshold (e.g. net per‑capita growth at low density with typical numbers of positive links). Even a rough analytical back‑of‑the‑envelope estimate would help readers see why a threshold should exist and how it depends on interaction benefits and self‑regulation.\n3. The contrast between strong vs weak inheritance is central to the eco‑evolutionary framing, but the current description leaves some questions.\n• The choice of Δ = 5 vs Δ = 30 feels somewhat ad hoc. While the qualitative distinction is clear, it would be helpful to comment on why these values were chosen and whether intermediate values (e.g. Δ = 10 or 15) behave similarly to one of these extremes.\n• Please clarify exactly how interaction changes are implemented when the parent has few links: is the number of deletions capped by the number of existing interactions? How is the balance of addition vs deletion chosen for a given d?\n• The model assumes inheritance acts only on interaction structure, not on intrinsic growth rates or self‑regulation. This is a reasonable simplification, but it is worth stating explicitly in the Discussion that “evolution” here acts purely on interactions, which may limit the analogy to full trait or demographic evolution.\n4. The Methods define an “effective increase” of degree entropy and modularity as (Z − Z_r)/Z_r, where Z_r is the mean for random networks. In the text, this is described as a “z‑score metric”, but it is not a z‑score in the usual statistical sense ((Z − mean)/SD). This should be renamed to avoid confusion.\nIn addition:\n• Both entropy and modularity are computed on the unweighted interaction network, seemingly ignoring sign and direction. Yet the underlying ecology is driven by signed, directed interactions (consumer vs resource, + vs –). Please state explicitly:\n– Whether the network is symmetrised and unsigned for these calculations.\n– Why this choice is appropriate for the questions at hand, and what information might be lost.\n• Since modularity in mutualistic or bipartite systems is often computed with more specialised methods, a brief acknowledgement of these simplifications would be useful, even if exploring signed or bipartite modularity is beyond the scope of this paper.\nMinor comments\n1. Fig. 1: The colour scales are clear, but please ensure that the palette is friendly to colour‑blind readers. Also, it might help to mention in the caption that each heatmap cell is based on 15 replicate simulations.\n2. Methods\n• Explicitly state how connectance is computed (e.g. C = L/S² on the binary directed network).\n• Mention that the extinction threshold (10⁻⁶) has been checked and does not qualitatively affect results, if that is the case.\n• As noted above, please rename the “z‑score” to something like “relative increase over random networks” to align with the formula.\nIn summary, this manuscript offers a valuable eco‑evolutionary perspective on how multiple interaction types, inheritance, and ecological filtering jointly shape community structure and macroecological patterns. Addressing the points above, especially clarifying the dynamical model, sharpening the interpretation of the T1/T2 threshold, and strengthening the macroecological comparison, would substantially enhance its clarity, robustness, and impact.", "reviewer_3": "This manuscript investigates the assembly of interaction networks, and how this process is shaped by inheritance and selection of interaction types. The study uses a Lotka-Volterra-type population dynamics model incorporating competitive, consumer-resource, and mutualistic interactions, and simulates scenarios in which new species are added either with or without inherited interactions. The main findings are: (1) the relative abundance of the three interaction types depends on the balance between interaction strength and interaction cost, (2) both inheritance and selection can promote network complexity, and (3) comparison with empirical microbial data suggests mechanisms that could have produced observed communities.\nThe novelty of this work lies in incorporating multiple interaction types, in contrast to most existing models that consider only one. This opens the door to addressing interesting questions in community assembly. The comparison with empirical data is also a strength. The manuscript is clearly written, the methodology is sound, and the figures are clear. However, several points would benefit from clarification, as I detail below.\nComment 1: overlap with previous study\nThere is substantial overlap with a previous publication (Araujo & Lurgi, 2025, “Mutualism provides a basis for biodiversity in eco-evolutionary community assembly”. PLoS Comput Biol.), which not acknowledged.\nFirst, both papers use the same eco-evolutionary model. While it is reasonable to apply the same model to explore different questions, the authors should acknowledge this explicitly (e.g., “We used a previously developed model…” in L119 rather than “We developed…”).\nSecond, the published study focused on scenarios with high strength-to-cost ratios that yielded predominantly mutualistic (T2) communities. The present manuscript expands the parameter space to explore conditions generating predominantly competitive (T1) or consumer-resource communities. However, the results concerning network complexity in mutualism-dominated communities (L184–189) and the role of evolution (L258–272) appear to be essentially the same in both studies. Agan, this is not acknowledged.\nIn summary, the manuscript should clearly outline similarities and differences with the previous study or reduce repetition where possible.\nComment 2: rationale for inheritance of interactions\nThe mechanism by which species “inherit interactions” from parents may seem abstract. One could view this as inheritance of traits that cause offspring to interact with similar partners as their parents. The mechanism itself is reasonable, but providing more biological justification or framing in the Introduction would help readers understand this modelling choice.\nComment 3: model details\nThe following clarifications would improve the model description:\n(a) It may be helpful to describe in words what each term in Equation 1 represents (e.g., “The first term represents intrinsic growth…” etc.).\n(b) My understanding is that interspecific interaction terms (competitive, mutualistic, consumer-resource) are nonzero only for interacting species pairs (i.e., cij, mij, pij = 0 for non-interacting pairs). Please clarify whether this is correct.\n(c) It is unclear how initial interactions are defined. My understanding is: at the start there are five species with no interactions; under invasion, a new species is introduced and is connected to some existing species according to the drawn connectivity; under evolution, a new species inherits interactions from an existing species. If this is the case, how can interactions be inherited if the initial five species have none?\n(d) My understanding is that there is no rewiring of interactions during ecological dynamics or after species addition. Species may lose interactions if their partners go extinct, which is the selection mechanism, but interactions are otherwise static. If correct, this simplification should be acknowledged, as rewiring is common in natural systems.\n(e) Including a conceptual figure of the model and simulation workflow would be helpful.\nComment 4: effect of intraspecific competition\nThe result in L189-192 and Figure S1 is interesting, but discussed only briefly. For example, T1 communities arise when sigma/delta is low and intraspecific competition (s) is low (as in the main text), but also when sigma/delta is high and s is high. Could the authors speculate on why the magnitude of intraspecific competition reverses the effect of the strength-to-cost ratio?\nComment 5: complexity in T1 communities\nThe discussion focuses on the emergence of complexity in T2 communities (L341–343, L369-370). Could the authors also comment on complexity in T1 communities? In contrast to T2 communities, the Evo (delta=5) and Inv.Prop. scenarios produce higher complexity than Evo (delta=30) and Inv.Rand.\nComment 6: comparison with empirical data\nThis comparison is an interesting addition. However, simulation results do not appear highly sensitive to the eight scenarios simulated. Have the authors examined sensitivity to other model parameters, such as intrinsic growth rate or intraspecific competition (which seems to have substantial effects elsewhere)?\nSpecific comments\n- L23: stability was not evaluated in this study and should be removed from this sentence\n- L275-276: the phrase “exhibit more structured properties” is somewhat vague, please consider clarifying what is meant"}, "author_response": "R. We thank the reviewer for the positive assessment of our manuscript and for the constructive suggestions. We have addressed all points raised and revised the manuscript accordingly. In particular, several of the reviewer’s suggestions motivated additional simulations and robustness analyses, substantially expanding the scope of scenarios explored and strengthening the generality of our conclusions.\nHarvesting cost delta: this seems a relevant parameter in your model and investigation, but it is not very well explained or justified. Why is it needed? Is there any empirical or theoretical foundation for it? Why is it constant? The name itself is also not very informative.\nR. We agree that the rationale behind this parameter required clearer explanation. We have expanded the description of the harvesting cost δ in the main text (lines [188–200]) to provide a stronger biological motivation:\n“[...] incurs a harvesting cost 𝛿, representing the energetic and physiological investment required to acquire and maintain benefits from a given interaction. Such costs arise from the allocation of finite resources to interaction-related traits. For simplicity and interpretability, we assume a constant per-interaction cost, understood as a mean-field approximation capturing the average energetic burden of maintaining a positive interaction. This cost is subtracted from 𝑟 for each beneficial interaction, reflecting a trade-off: as a species accumulates positive links, 𝑖 it becomes increasingly reliant on its partners and may ultimately depend on them for survival (i.e., obligate mutualism or increasing trophic dependence on external resources for consumers, with reduced viability in their absence). This choice reflects the allocation of energy in either interaction traits or survival and reproduction. The balance between the interaction strength parameter 𝜎 and the per‐link cost 𝛿 therefore determines the net benefit of positive interactions.”\nNew interactions for both evolutionary and invasion are extracted from a half-normal. The invasion scenarios Rand. and Prop. specify how the types of interactions are chosen, but it is unclear for the evolutionary scenarios.\nR. Thank you for pointing this out. We have clarified this in the manuscript in lines [261-263], as part of the definition of evolution scenarios:\n“The types of newly added interactions are determined by choosing the signs of each directed interaction (positive or negative) with equal probability.”\nIntroduction of new species happens at equilibrium, are results robust if mutations and/or invasions were not rare?\nR. We tested the robustness of our results to violations of the rare-event assumption by running additional simulations in which mutation (Evo 5) and invasion (Inv Prop) events occurred before ecological equilibrium was reached. These results are now presented in Supplementary Figures S9–S12. Overall, the qualitative patterns remained unchanged. Deviations were observed only in extremely fast assembly regimes, where limited relaxation time prevents communities from fully settling into their characteristic structural properties. We added this to our results in lines [336-337].\nEvolution and invasion are treated separately, what results would emerge if both were acting at the same time?\nR. We explored this by running additional simulations combining evolutionary (Evo 5) and invasion-based (Inv Prop) assembly events in varying proportions, now shown in Supplementary Figures S13 and S14. These simulations reveal gradual transitions between the patterns observed under pure evolution and pure invasion, with some nonlinear behaviour emerging when evolutionary events dominate. We highlight this in lines [336-342].\nAmong the network metrics, nestedness is not mentioned, but often used and discussed in this type of research.\nR. We agree that nestedness is often discussed in the context of mutualistic communities. However, nestedness is typically defined for bipartite networks. Because our communities involve multiple interaction types and are not naturally bipartite, analysing nestedness would require additional assumptions about how to partition interactions. For this reason, we chose not to include nestedness metrics in the present study.\nAlso, unweighted metrics are used, are results robust to the inclusion of weights in network metrics?\nR. We examined the robustness of our results to the inclusion of interaction weights, and present the corresponding analyses in Supplementary Figure S24. Most qualitative patterns remain unchanged. One notable difference is that weighted degree entropy in Type 1 communities shows a consistent decrease relative to random networks across all scenarios. Aside from this effect, the main conclusions are robust to the inclusion of weights. We report this in our results in lines [358-360].\nComparisons with empirical data are qualitative, but the similarity between simulations and data should rather be quantified formally.\nR. We agree that a formal quantitative comparison strengthens the empirical relevance of the study. We have therefore extended the data analysis to include quantitative comparisons between simulations and microbial data, as well as control simulations without interactions between species. Specifically, we computed Wasserstein distances between probability distributions of abundance correlations and present these results in the new supplementary section “Complementary data analysis scenarios (Figs. S25–S29)”. This analysis confirms that Type 2 communities are generally more consistent with empirical macroecological patterns than Type 1 communities, a conclusion now explicitly stated in the revised main text (lines [400–416]).\nI hope these suggestions can help improve the manuscript.\nR. We thank the reviewer again for these helpful suggestions, which significantly improved the clarity, scope, and robustness of the manuscript.\nReviewer #2 (Remarks to the Author):\nThe eco-evolutionary assembly of complex communities with multiple interaction types\nIn this study, the authors develop an eco-evolutionary community assembly model that integrates three interaction types (competition, consumer–resource, and mutualism) within a single dynamical framework. Communities assemble via sequential introduction of species, either through speciation with inheritance of interaction structure or via invasion from a regional pool. A cost–benefit trade‑off is imposed on each positive interaction (mutualistic or consumer) via a per-link harvesting cost and a saturating Type II functional response for the benefits. By exploring the plane of interaction strength (σ) and cost (δ), the authors identify a threshold separating two regimes: a competition-dominated, sparsely connected “Type 1” community and a mutualism-dominated, more connected “Type 2” community.\nThe paper then compares eight assembly scenarios that combine the two regimes (T1/T2), strong vs weak inheritance (Δ small vs large), and invasion with either fixed or variable interaction-type composition. This is used to separate the effects of selection of interaction types from inheritance of interactions on emergent complexity (richness × connectance), degree entropy, and modularity. Finally, the authors compare time series generated from assembled communities to empirical macroecological patterns from human microbiome datasets, focusing on mean abundance distributions and pairwise abundance correlations.\nUpon review, I am convinced that combining multiple interaction types with explicit inheritance in a community-assembly framework is timely and potentially impactful, especially given the connection to microbial macroecology. Generally, the model is well motivated and the results are interesting. However, I believe the manuscript can be substantially strengthened both conceptually and in clarity if the following major and minor points are addressed.\nR. We thank the reviewer for the careful and thoughtful evaluation of our manuscript. We are pleased that the reviewer finds the framework timely and the results interesting. We have addressed all major and minor comments and believe that the revisions substantially improve the conceptual clarity, transparency of the model, and interpretability of the results.\nMajor points\n1. The dynamical system is central to the paper, yet Eq. (1) is currently difficult to parse. As typeset, the different sums and terms (competition, consumer–resource, mutualism, costs) are condensed in a way that is hard to reconstruct, and some notation (e.g. which sums appear in the denominators, the meaning of “Σ_{p+,m} δ”, etc.) is not fully explained.\nI suggest:\n• Rewriting Eq. (1) in a more modular way (e.g. grouping explicitly into “self-regulation”, “loss to consumers”, “gain from resources”, “gain from mutualists”, “harvesting cost” terms), possibly split over several lines.\n• Explicitly defining all summation indices: over which sets (consumers, resources, mutualists) each sum runs, and which interactions contribute to the denominators of the Type II functional responses.\nR. We agree that the original presentation of Eq. (1) could be difficult to parse. We have therefore rewritten the equation in a more modular form and explicitly grouped terms according to their ecological roles (self-regulation, competitive effects, consumer-resource losses and gains, mutualistic gains, and intrinsic growth with costs). In addition, we now clearly define the scope of each summation and explicitly state which interaction terms contribute to the denominators of the Type II functional responses. These changes are described in the main text (lines [166–173]) and are further supported by a new conceptual figure (Fig. 1B), which summarises all mechanistic components of the dynamical equation.\n• Clarifying in the main text how consumer–resource asymmetry is handled (the rule p_{ij}⁺ ≤ p_{ji}⁻ and the enforced symmetry when this fails). This is an important modelling choice and deserves a brief justification, not only a note in Methods.\nR. We agree that this modelling choice deserves explicit justification in the main text. We have therefore expanded the explanation in lines [178–184]:\n“In the case of consumer-resource interactions, we impose a directional constraint on interaction strengths to reflect energetic inefficiency. Specifically, we assume that the number of resource individuals removed through consumption always exceeds the number of new consumer individuals produced as a consequence of that consumption. Accordingly, when interaction strengths are randomly sampled, if 𝑝+ > 𝑝−, we reset 𝑝+ so that 𝑝+ = 𝑝−. This 𝑖𝑗 𝑗𝑖 𝑖𝑗 𝑖𝑗 𝑗𝑖\nconstraint prevents unrealistically efficient trophic conversion.”\n• Stating explicitly whether the network is treated as directed and signed in the dynamics but undirected/unsigned when computing network metrics (see point 5).\nR. We now explicitly clarify this point in the manuscript (lines [347–349]): “[…] we analysed the topological degree entropy (a measure of the evenness in the number of links across species) and network modularity in assembled communities. For this analysis, community networks were represented in an undirected and unweighted form, retaining only the presence or absence of interactions”\n2. The authors show a sharp transition in interaction-type composition, connectance, and “complexity” as σ and δ vary, and then select a single parameter pair in each region as “Type 1” and “Type 2”. Supplementary Fig. S1 further demonstrates that rescaling intraspecific competition shifts communities across this threshold, indicating that it reflects the balance of interaction benefits vs self-regulation, rather than σ and δ alone.\nI suggest:\n• Emphasising that T1 and T2 are regimes or “community types” defined by this relative balance, not single parameter points. The selected (σ, δ) pairs should be clearly described as representative.\n• Providing at least a brief robustness check (in the SI) showing that key results for the assembly scenarios (e.g. the roles of inheritance and selection in Fig. 2–4) are qualitatively stable under modest changes in σ and δ within each regime.\nR. We agree that Type 1 and Type 2 communities should be interpreted as regimes defined by the balance between interaction benefits and costs, rather than as single parameter combinations. We have clarified this interpretation in the revised manuscript and explicitly describe the chosen (𝛿,𝜎) pairs as representative points within broader regions of parameter space (lines [284–289]):\n“Community types (T1 and T2): for every scenario, we chose particular points (δ,σ) to represent T1 and T2 community regions, (δ = 2.5∗10−2,σ = 0.05) for T1 (fifth column from the left and second row from the bottom in Fig 2) and (δ = 1.0∗10−2,σ = 0.2) for T2 (second column from the left and fifth row from the bottom in Fig 2). We tested different points within the same regions and all yielded similar analyses (Figs S5-S8).”\n• Considering a simple heuristic or mean-field argument for the location of the threshold (e.g. net per-capita growth at low density with typical numbers of positive links). Even a rough analytical back-of-the-envelope estimate would help readers see why a threshold should exist and how it depends on interaction benefits and self-regulation.\nR. We thank the reviewer for this suggestion. We have added a simple mean-field analysis that considers the net per-capita gain of a focal species acquiring an additional positive interaction at equilibrium density. This analysis reveals a threshold in the space defined by the number of positive interactions and the per-link cost, beyond which accumulating additional positive interactions becomes detrimental. This threshold depends on self-regulation, which controls equilibrium abundances and therefore modulates the benefits of positive interactions. This analytical argument, presented in the main text (lines [224–243]) and Supplementary Figure S1, provides intuition for the emergence of the two community regimes observed in simulations.\n3. The contrast between strong vs weak inheritance is central to the eco-evolutionary framing, but the current description leaves some questions.\n• The choice of Δ = 5 vs Δ = 30 feels somewhat ad hoc. While the qualitative distinction is clear, it would be helpful to comment on why these values were chosen and whether intermediate values (e.g. Δ = 10 or 15) behave similarly to one of these extremes.\nR. We agree that the choice of Δ values requires justification. We therefore extended the analysis to include intermediate inheritance strengths (Δ = 14 and 22), now presented in Supplementary Figures S3 and S4. These intermediate cases behave qualitatively closer to the weak-inheritance regime (Δ = 30), indicating that our conclusions are not sensitive to the precise numerical values chosen for Δ in this case and that Δ=5 successfully captures a strong-inheritance case. We present these findings in lines [266-267].\n• Please clarify exactly how interaction changes are implemented when the parent has few links: is the number of deletions capped by the number of existing interactions? How is the balance of addition vs deletion chosen for a given d?\nR. We now clarified this in lines [136-140]: “For a given 𝑑, the number of interactions created is uniformly chosen from [1,𝑑] and the number of interactions destroyed is 𝑑 minus that number. If there are not enough interactions inherited (or no interactions, as in the initial state), then new interactions are created whenever an inexistent one would be destroyed (and vice-versa).”\n• The model assumes inheritance acts only on interaction structure, not on intrinsic growth rates or self-regulation. This is a reasonable simplification, but it is worth stating explicitly in the Discussion that “evolution” here acts purely on interactions, which may limit the analogy to full trait or demographic evolution.\nR. We agree that this assumption should be stated explicitly. We now clarify in the Discussion that evolution is about interactions both in lines [427], “Evolution of interactions by speciation […]”, and [483-484], “[…] evolutionary processes (acting on species interactions) […]”.\n4. The Methods define an “effective increase” of degree entropy and modularity as (Z − Z_r)/Z_r, where Z_r is the mean for random networks. In the text, this is described as a “z-score metric”, but it is not a z-score in the usual statistical sense ((Z − mean)/SD). This should be renamed to avoid confusion.\nR. We agree and have removed the labelling of a z-score metric, now stating “[…] a score metric of relative increase over random networks […]”, in lines [352-353].\nIn addition:\n• Both entropy and modularity are computed on the unweighted interaction network, seemingly ignoring sign and direction. Yet the underlying ecology is driven by signed, directed interactions (consumer vs resource, + vs –). Please state explicitly:\n– Whether the network is symmetrised and unsigned for these calculations.\n– Why this choice is appropriate for the questions at hand, and what information might be lost.\n• Since modularity in mutualistic or bipartite systems is often computed with more specialised methods, a brief acknowledgement of these simplifications would be useful, even if exploring signed or bipartite modularity is beyond the scope of this paper.\nR. We agree that this choice required clearer justification. We now explicitly state that degree entropy and modularity are computed on symmetrised, unsigned, and unweighted versions of the interaction networks. We clarify that this choice is deliberate, as our aim is to characterise the topological interaction architecture. We also explicitly acknowledge the information discarded by this simplification. These clarifications have been added in lines [343–351]:\n“To further examine how interaction architecture (i.e. the structure of who interacts with whom) emerges during community assembly, independently of interaction sign, direction, or magnitude, we analysed the topological degree entropy (a measure of the evenness in the number of links across species) and network modularity of assembled communities. For this analysis, community networks were represented in an undirected and unweighted form, retaining only the presence or absence of interactions. This choice isolates structural properties of interaction architecture that are robust to dynamical fluctuations in interaction strengths, while necessarily discarding information on interaction type, directionality, and energetic asymmetries.”\nMinor comments\n1. Fig. 1: The colour scales are clear, but please ensure that the palette is friendly to colour-blind readers. Also, it might help to mention in the caption that each heatmap cell is based on 15 replicate simulations.\nR. We confirm that all colour palettes used in the figures were selected to be colourblind-friendly. We have ensured consistency across figures. The figure caption (now Figure 2) specifies that the results presented are for 15 replicated simulations.\n2. Methods\n• Explicitly state how connectance is computed (e.g. C = L/S² on the binary directed network).\nR. We now state it as suggested.\n• Mention that the extinction threshold (10⁻⁶) has been checked and does not qualitatively affect results, if that is the case.\nR. We have tested and now explicitly stated that our particular choice of threshold does not qualitatively affect results, in lines [557-560]:\n“A species is considered extinct and removed from the system if its abundance falls below a threshold of 𝑥 = 10−6 (we verified that results are qualitatively equivalent for 𝑥 = 10−5 𝑒𝑥𝑡 𝑒𝑥𝑡\nand 𝑥 = 10−7 and thus not sensitive to a particular choice).”\n𝑒𝑥𝑡\n• As noted above, please rename the “z-score” to something like “relative increase over random networks” to align with the formula.\nR. We renamed it as mentioned above.\nIn summary, this manuscript offers a valuable eco-evolutionary perspective on how multiple interaction types, inheritance, and ecological filtering jointly shape community structure and macroecological patterns. Addressing the points above, especially clarifying the dynamical model, sharpening the interpretation of the T1/T2 threshold, and strengthening the macroecological comparison, would substantially enhance its clarity, robustness, and impact.\nR. We thank the reviewer for the insightful and constructive comments. We believe that the revisions have substantially improved the clarity, interpretability, and conceptual strength of the manuscript.\nReviewer #3 (Remarks to the Author):\nThis manuscript investigates the assembly of interaction networks, and how this process is shaped by inheritance and selection of interaction types. The study uses a Lotka-Volterra-type population dynamics model incorporating competitive, consumer-resource, and mutualistic interactions, and simulates scenarios in which new species are added either with or without inherited interactions. The main findings are: (1) the relative abundance of the three interaction types depends on the balance between interaction strength and interaction cost, (2) both inheritance and selection can promote network complexity, and (3) comparison with empirical microbial data suggests mechanisms that could have produced observed communities.\nThe novelty of this work lies in incorporating multiple interaction types, in contrast to most existing models that consider only one. This opens the door to addressing interesting questions in community assembly. The comparison with empirical data is also a strength. The manuscript is clearly written, the methodology is sound, and the figures are clear. However, several points would benefit from clarification, as I detail below.\nR. We thank the reviewer for the careful, thorough, and constructive evaluation of our manuscript. We particularly appreciate the attention given to issues of transparency, biological interpretation, and model clarity.\nComment 1: overlap with previous study\nThere is substantial overlap with a previous publication (Araujo & Lurgi, 2025, “Mutualism provides a basis for biodiversity in eco-evolutionary community assembly”. PLoS Comput Biol.), which not acknowledged.\nFirst, both papers use the same eco-evolutionary model. While it is reasonable to apply the same model to explore different questions, the authors should acknowledge this explicitly (e.g., “We used a previously developed model…” in L119 rather than “We developed…”).\nSecond, the published study focused on scenarios with high strength-to-cost ratios that yielded predominantly mutualistic (T2) communities. The present manuscript expands the parameter space to explore conditions generating predominantly competitive (T1) or consumer-resource communities. However, the results concerning network complexity in mutualism-dominated communities (L184–189) and the role of evolution (L258–272) appear to be essentially the same in both studies. Agan, this is not acknowledged.\nIn summary, the manuscript should clearly outline similarities and differences with the previous study or reduce repetition where possible.\nR. We agree that this issue is important and thank the reviewer for raising it. The present manuscript builds on an eco-evolutionary modelling framework that we previously developed, and we now explicitly acknowledge this throughout the revised text. In particular, we have replaced phrasing suggesting model development with language clarifying that we use a previously developed model, and we cite the earlier study wherever appropriate. While both studies share the same modelling backbone, they address distinct research questions. The earlier paper focused on whether mutualism, under evolutionary assembly, promotes biodiversity, complexity, and stability, with an emphasis on mutualism-dominated communities. In contrast, the present manuscript investigates how the balance between interaction benefits and costs, together with inheritance and ecological selection of interaction types, gives rise to qualitatively distinct community regimes and determines their structural and macroecological properties. Some results for mutualism-dominated communities partially overlap between the two studies. However, in the present manuscript these results are embedded in a different conceptual framework and are used to address different questions, particularly regarding regime emergence, the separation of inheritance from ecological filtering, and the consequences for community architecture. We have revised the manuscript to ensure that these connections are explicit and that any overlap is appropriately acknowledged.\nComment 2: rationale for inheritance of interactions\nThe mechanism by which species “inherit interactions” from parents may seem abstract. One could view this as inheritance of traits that cause offspring to interact with similar partners as their parents. The mechanism itself is reasonable, but providing more biological justification or framing in the Introduction would help readers understand this modelling choice.\nR. We agree that this mechanism benefits from clearer biological framing. We included a better explanation as part of our definition of assembly by evolutionary speciation, before starting to present our results, in lines [126-129]:\n“Inheritance with changes in species interactions represents mutations in foraging, morphological, or physiological traits that biologically determine interactions with other species (e.g. feeding modes, body morphology, metabolic pathways, signaling mechanisms, etc.).”\nComment 3: model details\nThe following clarifications would improve the model description:\n(a) It may be helpful to describe in words what each term in Equation 1 represents (e.g., “The first term represents intrinsic growth…” etc.).\nR. We agree, and we have now added a detailed verbal description of each term in Eq. (1) in the main text (lines [166–173]). In addition, we introduced a new conceptual figure (Fig. 1B) that summarises all mechanistic contributions to the ecological dynamics.\n(b) My understanding is that interspecific interaction terms (competitive, mutualistic, consumer-resource) are nonzero only for interacting species pairs (i.e., cij, mij, pij = 0 for non-interacting pairs). Please clarify whether this is correct.\nR. Yes, this is correct. We clarified this in lines [160-163]: “The non-negative coefficients 𝑝 , 𝑖𝑗 𝑚 , and 𝑐 represent the strengths of consumer-resource (+/-), mutualistic (+/+), and 𝑖𝑗 𝑖𝑗 competitive (-/-) interactions, respectively, and are nonzero only when there is an interaction of the corresponding type between 𝑖 and 𝑗”.\n(c) It is unclear how initial interactions are defined. My understanding is: at the start there are five species with no interactions; under invasion, a new species is introduced and is connected to some existing species according to the drawn connectivity; under evolution, a new species inherits interactions from an existing species. If this is the case, how can interactions be inherited if the initial five species have none?\nR. This point was also mentioned by reviewer #2, and this is our response as stated above: we now clarified this in lines [136-140]: “For a given 𝑑, the number of interactions created is uniformly chosen from [1,𝑑] and the number of interactions destroyed is 𝑑 minus that number. If there are not enough interactions inherited (or no interactions, as in the initial state), then new interactions are created whenever an inexistent one would be destroyed (and vice-versa).”\n(d) My understanding is that there is no rewiring of interactions during ecological dynamics or after species addition. Species may lose interactions if their partners go extinct, which is the selection mechanism, but interactions are otherwise static. If correct, this simplification should be acknowledged, as rewiring is common in natural systems.\nR. The reviewer is correct in this. We now acknowledge this simplification in lines [155-157]:\n“After a new species is introduced in an assembly event, we assume no rewiring of its interactions, only the possible addition of interactions with new species included in subsequent events.”\n(e) Including a conceptual figure of the model and simulation workflow would be helpful.\nR. This is a very helpful suggestion and have added a new conceptual figure (Fig. 1) illustrating the community assembly process and the mechanistic components of the model. Figure 1A depicts the sequence of assembly events and ecological equilibria, while Figure 1B summarises the contributions of each term in the dynamical equation.\nComment 4: effect of intraspecific competition\nThe result in L189-192 and Figure S1 is interesting, but discussed only briefly. For example, T1 communities arise when sigma/delta is low and intraspecific competition (s) is low (as in the main text), but also when sigma/delta is high and s is high. Could the authors speculate on why the magnitude of intraspecific competition reverses the effect of the strength-to-cost ratio?\nR. The reviewer is correct in their interpretation / assessment of the sigma/delta effect, however a higher intraspecific competition favours T1 communities, while a lower intraspecific competition favours T2 communities. We now made our statement clearer on this aspect by adding in line [217-220]:\n“we found that the intensity of intraspecific competition shifts the position of the threshold (Supplementary Fig S2), with higher values shifting the regime from T2 to T1, indicating that the benefits of interactions must be evaluated relative to the strength of self-regulation”.\nWe now also included a mean field exploration of the threshold (lines [224-243]). Our results indicate that the intraspecific competition changes the threshold by changing the mean species abundances, which in turn change how beneficial are the positive interactions.\nComment 5: complexity in T1 communities\nThe discussion focuses on the emergence of complexity in T2 communities (L341–343, L369-370). Could the authors also comment on complexity in T1 communities? In contrast to T2 communities, the Evo (delta=5) and Inv.Prop. scenarios produce higher complexity than Evo (delta=30) and Inv.Rand.\nR. We have expanded the Discussion, now addressing complexity in Type 1 communities, in lines [471-480]:\n“Moreover, we also highlight the role of the ecological selection of interaction types as a generator of complexity, since it enhanced the emergence of complexity in Type 1 communities in the scenarios with effective selection (evolution with high inheritance and invasion with variable proportion of interaction types). This happened regardless of these scenarios featuring the highest proportions of competition, which means that complexity emerges through selection not only with a joint domination of mutualism. Both species richness and connectance were on average higher for these scenarios (Fig 4B), although in Type 2 communities the effect of ecological selection of interaction types on complexity resulted from an increase in species richness only (Fig 4D, difference between the two invasion scenarios).”\nComment 6: comparison with empirical data\nThis comparison is an interesting addition. However, simulation results do not appear highly sensitive to the eight scenarios simulated. Have the authors examined sensitivity to other model parameters, such as intrinsic growth rate or intraspecific competition (which seems to have substantial effects elsewhere)?\nR. We added a clarification in line [383] that results are not sensitive to scenarios in T1 communities due to the sparsity of interactions. We ran extra simulations varying the levels of intrinsic growth rate and intraspecific competition and presented results in Supplementary Figures S15-S22, and also mixed evolution and invasion scenarios presented in Figures S13 and S14. For the comparison with data, we presented the most interesting results among all these extra simulations in Figures S28 and S29, all for T2 communities. In addition to including some extra scenarios, we quantified the comparison with data, for the correlation pattern, using the Wasserstein distance between probability distributions, comparing data, simulations, and the logistic counterparts with no interactions between species. In this way, we could better assess the specific contribution of interactions in making the simulated pattern more similar to the data. The additional analyses made more robust our conclusion that T2 communities result in a better match to empirical patterns. However, it is hard to point out a better scenario overall, since there are many aspects to consider and the patterns were indeed similar across all scenarios, even considering the extra simulations. All these new analyses and quantitative analyses of empirical data matching are presented in lines [336-342] and lines [400-416].\nSpecific comments\n- L23: stability was not evaluated in this study and should be removed from this sentence\nR. We agree and have removed this statement.\n- L275-276: the phrase “exhibit more structured properties” is somewhat vague, please consider clarifying what is meant\nR. We agree that this phrasing was vague and have removed it."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The Authors carefully addressed all my concerns, and as far as I can tell also those from other Referees. I do not have any additional major considerations, however I attach a few minor suggestions on the main text and supplementary materials.", "reviewer_2": "I am happy with the revised manuscript. The authors have carefully addressed all of my comments, and the revisions have clearly improved the clarity and presentation of the model, as well as the interpretation of the results.\nOverall, the manuscript is now clear, well motivated, and ready for publication.", "reviewer_3": "I thank the Authors for addressing my comments and answering my questions. I do not have any further comments. I believe this work to be an interesting contribution to the field."}, "author_response": "R. We thank the reviewer for their careful evaluation of the revised manuscript and for confirming that the main concerns have been satisfactorily addressed. We appreciate their additional minor suggestions and have incorporated them to the revised version of the manuscript. We thus have improved the conceptual figure, the presentation of equations, and included the Wasserstein distances in Fig 5.\nR. We thank the reviewer for their constructive feedback throughout the review process. We are pleased that the revisions have improved the clarity of the model and the interpretation of the results, and we appreciate the reviewer’s positive evaluation of the final version.\nR. We thank the reviewer for their thoughtful comments and questions, which helped us strengthen the manuscript. We appreciate their positive assessment and their recognition of the contribution of this work to the field."}]

Decision letter