The joint role of coevolutionary selection and network structure in shaping trait matching in mutualisms

Coevolution can sculpt remarkable trait similarity between mutualistic partners. Yet, it remains unclear which network topologies and selection regimes enhance trait matching. To address this, we simulate coevolution in topologically distinct networks under a gradient of mutualistic selection strength. We describe three main insights. First, trait matching is jointly influenced by the strength of mutualistic selection and the structural properties of the network where coevolution is unfolding. Second, the strength of mutualistic selection determines the network descriptors better correlated with higher trait matching. While network modularity enhances trait matching when coevolution is weak, network connectance does so when coevolution is strong. Third, the structural properties of networks outrank those of modules or species in determining the degree of trait matching. Our findings suggest networks can both enhance or constrain trait matching, depending on the strength of mutualistic selection.


Introduction
All organisms are subjected to biological evolution. Yet, it is argued, much of evolution is in fact coevolution. That is, the reciprocal evolutionary adaptation between interacting species mediated by natural selection [1]. From shaping remarkable trait complementarity between interacting partners [2][3][4] to driving the major evolutionary transitions [5], coevolution is an organizer of biodiversity.
While mainly documented in pairwise interactions [6], coevolution is not restricted to this scale. Rather, it operates on ecological networks; though the extent to which coevolution sculpts networks is contended. While complex network structures have been found to arise from coevolutionary processes [7][8][9][10], they can also emerge without a shared evolutionary history [11,12]. Hence, the recurrent topologies of networks found in nature [13,14] likely result from a suite of different mechanisms including a potential role of coevolution as shown experimentally with bacteria and phages [9].
Beyond being potentially organized by coevolution, the topology of ecological networks can in turn influence coevolutionary dynamics [9,15,16]. The wiring of networks impacts how selective pressures ripple across communities and ultimately affects trait patterns [8,[16][17][18]. A major challenge lies in understanding and predicting the dynamics and traits generated by coevolution when scaling up from the level of pairwise interactions to complex networks [19]. In this respect, initial theoretical studies have shown indirect effects can be as important as direct interactions in driving the coevolution of traits in large mutualistic networks. [16]. Furthermore, Guimarães et al. [16] found the contribution of indirect effects to coevolution differed across structures of interaction networks. Thus, network architecture seems to be an important regulator of coevolutionary dynamics.
Much emphasis has been placed on the interplay between network structure, ecological [20][21][22], and evolutionary processes [16,23]. In the case of coevolution, selection imposed by interacting partners (i.e. coevolutionary selection) is required. Thus, the strength of coevolutionary selection will regulate the evolutionary consequences of interactions. In fact, in mutualistic networks, the strength of coevolutionary selection (i.e. mutualistic selection) can determine the degree of trait matching between interacting partners [8,18,24] and how redundant interactions are in a network [18]. Unravelling how network properties and the strength of mutualistic selection jointly influence the emergence of trait similarity between partners is particularly relevant as higher trait matching and redundancy have been shown to increase the stability of mutualistic networks [18].
Here, we build upon a model of coevolution in networks [16] to explore how network architecture and the strength of mutualistic selection affects the evolution of trait matching between partners (figure 1). We simulate coevolution in a set of topologically distinct artificial networks under a gradient of mutualistic selection strength and measure the degree of trait matching between partners after evolution. In doing so, we aim to answer the question: which network topologies and selection regimes result in greater trait matching?

Methods (a) Generating networks
We constructed all binary bipartite networks using a simple niche-based consumer-resource model [25]. We chose this model for two reasons. First, it has been shown to accurately reproduce the structural characteristics of mutualistic networks [25]. Second, networks are built with a pre-specified size and a target connectance. Thus, using this model, we can generate a set of comparable networks that encapsulate topological diversity but maintain the structural characteristics of real-life mutualistic networks.
The algorithm for building networks is based on the niche model initially proposed by Williams & Martinez [26], but modified to allow the building of bipartite networks [25]. In short, the model classifies species into two types: consumers and resources. It then establishes interactions between species if the niche of a consumer overlaps with the niche of a resource. The algorithm for building networks consists of four steps. First, all species in the community are initially assigned a random niche value n i from a uniform distribution on the interval [0, 1]. Second, a range value r i is defined for all consumer species, where r i = xn i and x is a random variable drawn from a beta-distributed probability function defined as p(x) = β(1 − x)(β − 1) and β = (1/2c). Third, a range centre c i is defined as a random value from a uniform distribution of range r i /2 and min(n i , 1 − r i /2). Fourth, interactions are established between a consumer species (i) and all of the resource species that fall into its diet interval We used the niche network model to build networks consisting of 40, 60, 80, and 100 species. For each network size, we constructed networks with different consumer and resource species ratios (table 1). For each ratio, we built networks with target connectance that ranged from 0.05 to 0.49 with increments of 0.02. For each value of connectance, we built 20 networks. In total, we constructed 5520 networks.

(b) Network, module, and species descriptors
We measured the modularity of each network using the Q metric [27], which relies on the identification of modules.  Figure 1. The interplay between network structure and mutualistic selection in driving trait evolution. We express communities in terms of their interaction networks. Here, we show two representative examples of communities that differ in their network structure (compare network structure of communities A and B). We simulate coevolution in the communities and record the trait values of all species in the network. We measure the degree of trait matching, represented here as a heatmap relating trait values of interacting species across guilds, at the network, module, and species level. We explore how differences in the level of trait matching arise due to the structural dissimilarities of our communities (note differences in colouration of the heatmap between communities) and due to the strength of mutualistic selection driving coevolution (compare heatmaps of the same community for the three levels of mutualistic selection). (Online version in colour.) royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211291 sets of highly connected species that are loosely connected to species from other modules. For a given network partition into modules, Q is computed as the difference between the observed fraction of interactions between species in the same module and the expected fraction of interactions connecting species in the same module if interactions were established at random [27,28]. Thus, Q = 0 if the number of observed interactions in modules is equally distributed compared with the expected number of interactions in modules when links are random. We used the MODULAR software [28] with a simulated annealing algorithm [29] to find the network partition where Q is maximized. For each network, we recorded the network partition that maximized Q and the module membership of each species. For each network partition, we measured seven topological properties. These descriptors capture the structural diversity at the network, module, and species scale. At the network scale we recorded: (i) network modularity (Q), (ii) network connectance (N c ), and (iii) the number of modules in the network (N m ). At the scale of modules, we recorded: (iv) module size (M s ) and (v) module connectance (M c ). Finally, at the species scale, we measured two properties that summarize their role in their respective networks: (vi) within-module degree (WMD) and (vii) among-module connectivity (AMC) [29,30]. The standardized WMD (zd) measures how well a species is connected to other species in its module. It is defined as follows: where k i s is the number of links species i has to other species in its module s, k s and σ ks are the mean and standard deviation of the within-module links of all species in s. The AMC c i measures how a species i connections are distributed between its own and other modules. It is defined as follows: where k i is the number of links species i has, k i s is k i in module s, and N M is the number of modules in a network. c i is bounded from 0 to 1. When c i = 0 all links of species i are contained in its module; when c i = 1 the links of species i are evenly distributed between all modules. The networks we built had different structures, yet shared similar topologies with empirical mutualistic networks (electronic supplementary material, figure S1). Q ranged from 0.08 to 0.81 with a mean of 0.37. N c ranged from 0.02 to 0.35 with a mean of 0.14. Networks had an average of 2.6 modules with a minimum of 1 and a maximum of 8. The smallest module was composed of two species, while the largest of 60. On average, modules contained 25 species. M c was greater than N c . M c ranged from 0.15 to 1, and on average modules had a connectance of 0.70. Following the classification of species roles by Olesen et al. [30], from our species pool (n = 129, 048, 840), 56% were peripherals (i.e. species with few interactions, mostly concentrated in their own module), 40.1% were connectors (i.e. species whose links connect different modules), 3% were module hubs (i.e. highly connected species with links predominately distributed inside their module), and 0.7% were network hubs (i.e. species acting as both connectors and module hubs).

(c) Coevolutionary model
We simulated coevolution using a discrete model that uses a selection gradient approach to link trait evolution with the fitness consequences of interactions [16]. Given a network of size S, the model takes N populations of species [1..S] and assigns a single trait value to each individual of each species. This trait value is assumed to affect both the fitness derived benefits from mutualistic interactions and from the environment. Thus, the mean trait value of the local population of species i (Z i ) evolves as a result of distinct mutualistic and environmental selective pressures. The evolution in discrete time of Z i is defined by where F is proportional to the selection gradient and additive genetic variance, N is the number of species in the network, q ij is the evolutionary effect of species j on species i, where 0 ≤ q ij ≤ 1, and θ i is the environmental optima for species i. Under this framework, selection due to interactions is assumed to favour trait matching between interacting species. The evolutionary effects of species j on species i are defined by ij ), α is a constant that determines the sensitivity of the evolutionary effect to trait matching, and a ij is an element of an adjacency matrix A of a given mutualistic network. a ij = 1 if species i interacts with species j, and a ij = 0 if they do not. The level of mutualistic selection parameter (0 ≤ m i ≤ 1) determines the relative importance of mutualistic selection in driving trait evolution. When m = 0, trait evolution is exclusively driven by environmental selection. When m = 1, trait evolution is exclusively driven by interactions.

(d) Simulating coevolution
To analyse how network structure affects trait evolution, we performed numerical simulations of the coevolution model (equations (2.3) and (2.4)) for all generated networks. Furthermore, in order to assess how the strength of mutualistic selection influenced the outcomes of trait evolution, we ran simulations for a gradient of m, ranging from 0.1 to 0.95 with increments of 0.05. We ran 20 replicas for each network at every level of mutualistic selection for a total of 110 400 simulations.
The values for θ, F, m, and the initial trait values for all species (z 0 i ) were sampled from a statistical distribution at the start of each simulation. We ran simulations until equilibrium, defined as j Z t i À Z tþ1 i j e. e = 1e − 6 and α = 0.02 as in Guimarães et al. [16]. The parameter values used in the simulations are described in electronic supplementary material, table S1. Though we drew trait values from U(0, 10), the results were qualitatively similar to those obtained when drawing them from N(5, 0.1), β(2, 5), or β(5, 2) (electronic supplementary material, figure S2). Furthermore, we set F to 0.5, a value that is greater than that used in Guimarães et al. [16]. F controls the slope of the selection gradient. Larger values increase the rate Table 1. Networks used to simulate coevolution. For each size, we built networks with different consumer (C) and resource (R) species ratios. For each ratio, we built networks with target connectance that ranged from 0.05 to 0.49 with increments of 0.02.

species
60 species 80 species 100 species royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211291 of adaptive change (electronic supplementary material, figure S3) and reduce the number of generations needed to reach equilibrium (electronic supplementary material, figure S4). By setting F ¼ 0:5 we significantly reduced the time required to perform the simulations, while not qualitatively affecting the overall results (electronic supplementary material, figure S5).

(e) Trait matching and null scenario
After each simulation, we calculated the level of trait matching between all species pairs belonging to different groups (i.e. consumers and resources) in the network. We defined trait matching (T) between a pair of species i, j at a given time t as For each network and at each level of m, we calculated the mean trait matching of each species (T s ), each module (T m ), and of the network (T n ) across replicas. We standardized the level of trait matching at the species (zTs), module (zTm), and network level (zTn) at each level of mutualistic selection. The standard trait matching score for a species i at a mutualistic level μ, was calculated as where T iμ is the trait matching level of species i at the mutualistic selection strength μ. T im and s Tm are the mean and standard deviation of trait matching for all species at μ. To contrast our findings from the coevolution model (hereafter coevolution scenario), with those expected by chance, we used a null coevolution scenario as in de Anderazzi et al. [10]. In this scenario, we fixed species traits to their environmental optima. We then calculated trait matching in the same way as in the coevolution model. We used the null coevolution scenario as a baseline for how similar species traits can be by chance, without coevolution operating.

(f ) Multivariate analyses
We used principal component analyses (PCA) to summarize structural properties into a single descriptor for each scale: network, module, and species. Hence, we fitted three PCAs using the relevant descriptors for each scale (see Network, module, and species descriptors section). We next describe the three PCAs. At the network scale, the first component of the PCA (hereafter referred to as network structure) explained 87% of variance. The network structure axis was positively associated with network modularity (0.96) and the number of modules in the network (0.89) and negatively associated with network connectance (− 0.95; electronic supplementary material, figure S6A). At the module scale, the first component of the PCA (hereafter referred to as module structure) explained 66% of variance. The module structure axis was positively associated with module size (−0.82) and negatively associated with module connectance (0.82, figure  S6b). At the species scale, the first component of the PCA (hereafter referred to as species role) explained 65% of variance. The species role axis was positively associated with both WMD (−0.81) and AMC (0.81, figure S6c). We placed each network, module, and species along their corresponding structural axis. We used a Spearman-Rank correlation to asses the relationship between network, module, and species descriptors with trait matching. We performed a correlation test at each level of mutualistic selection.
We used a redundancy analysis (RDA) to explore how the seven structural descriptors were related to species' trait matching. RDA is a type of constrained ordination analysis used to determine how much of the variation of a set of variables (response variables) is explained by the variation of another set of variables (explanatory variable, [31]). The RDAs we fitted used network connectance, number of modules, modularity, module size, module connectance, species' AMC, and WMD as response variables and species' standardized trait matching as the explanatory variable. We fitted an RDA at each level of mutualistic selection. The scores of the response variables on the RDA and on the first component of the PCA are shown in electronic supplementary material, tables S2 and S3.

Results
We first describe how network, module, and species structural descriptors individually influence the degree of trait matching arising at their corresponding scale. We contrast results of the coevolutionary model with those of the null scenario, where trait values were fixed to their environmental optima.
At the network scale, we found a positive correlation (ρ > 0) between network structure and trait matching for m ≤ 0.6, and a negative correlation (ρ < 0) for m > 0.6 (figure 2a). In other words, when mutualistic selection is weak (m < 0.6) the more modular a network is, the higher trait matching it attains. Yet, when mutualistic selection is strong (m > 0.6), network connectance is correlated with high trait matching. All correlations were statistically meaningful (α = 0.05) across values of mutualistic selection (electronic supplementary material, figure S7). Conversely, correlations were statistically meaningful in only one third of the levels of mutualistic selection levels for the null scenario (electronic supplementary material, figure  S8). Furthermore, ρ-values were much smaller (ranging from −0.08 to 0.07) compared to those of the coevolution scenario (ranging from −0.89 to 0.58, figure 2a). At the module scale, we found a negative correlation between module structure and trait matching across all levels of mutualistic selection (figure 2b). Note that the correlations were not statistically meaningful when 0.5 ≤ m ≤ 0.7 (electronic supplementary material, figure S9). By contrast, all correlations were statistically meaningful for the null scenario (electronic supplementary material, figure S10). Similarly, we found comparable ρ-values for both scenarios (ranging from −0.06 to − 0.05 for the null scenario and ranging from −0.11 to −0.0005 for the coevolution scenario; figure 2b). Hence, the apparent trend of large poorly connected modules attaining the highest levels of trait matching could be simply attributed to the structure of the networks and not to coevolution.
At the species scale, we found a negative correlation between role and trait matching when m ≤ 0.6 and a positive relation for m > 0.6 (figure 2c). Following Olesen et al. [30], a large WMD and AMC is characteristic of network hubs. Peripherals show the opposite trend. Thus, when mutualistic selection is weak (m < 0.6), peripherals attain the highest levels of trait matching. When selection is strong (m > 0.6), network hubs record the highest levels of trait matching. All correlations were statistically meaningful (α = 0.05) across mutualistic selection (electronic supplementary material, figure S11). While correlations were statistically meaningful in 44% of the mutualistic selection levels for the null scenario (electronic supplementary material, figure S12), ρ-values were much smaller (ranging from −0.004 to 0.004) than those of the coevolution scenario (ranging from −0.11 to 0.27; figure 2c).
We next describe how network, module, and species structural descriptors jointly influence the degree of trait matching arising at the scale of species. In other words, we explore how the level of trait matching a species achieves relates to the role it plays in its network, the structure of the module it belongs to, and the topology of the network it belongs to. We found trait matching explains more variation in structural descriptors as mutualistic selection increases (figure 3a). These patterns do not emerge for the null coevolution scenario. Furthermore, we found network modularity and connectance to be the predominant structural descriptors shaping trait matching ( figure 3b). Yet, the relation between network structure and trait matching changes with the strength of mutualistic selection. Higher trait matching is associated with increasing network modularity up to m ≤ 0.5 (figure 3b). When m > 0.65, higher trait matching is associated with increasing network connectance (figure 3b).

Discussion
Our results highlight three main insights. First, the degree of trait matching arising as a result of coevolution is jointly influenced by the strength of mutualistic selection and the structural properties of the network where coevolution is unfolding.
Second, the set of network descriptors correlated with higher trait matching is determined by the strength of mutualistic selection. Third, the structural properties of networks outrank those of modules or species in determining the evolved degree of trait matching. We expand upon these points below.
Stronger mutualistic selection results in greater trait matching (electronic supplementary material, figures S13-S15). This is expected. If mutualistic selection increases trait complementarity between partners, and the strength of selection is intensified, this will result in greater trait matching [8,[16][17][18]. We found this trend to occur regardless of the scale at which trait matching is measured (either at the scale of networks, modules, or species). However, because mutualistic selection traverses networks, their structure will play a role in how coevolution operates [16]. We found network structure influences the degree of trait matching evolved for most levels of mutualistic selection. Strikingly, however, when mutualistic selection is equally important as environmental selection in driving trait evolution, we found no association between structure and degree of trait matching. Thus, for the most part, the observed coevolved trait complementarity in partners is a joint effect of how strong mutualistic selection is and what the structural properties of a network are.
The network descriptors associated with greater trait matching change with the strength of mutualistic selection. At the network scale, we do not find modularity or connectance to consistently increase trait matching in a community. Both amplify trait complementarity between partners but under different circumstances. When mutualistic selection is weak, network modularity enhances trait matching in a community. By contrast, network connectance amplifies trait matching when mutualistic selection is strong. These differences highlight how network architecture can modulate the role of mutualistic selection in shaping trait evolution. When mutualistic selection is weak, high trait complementarity between partners will likely only arise within tightly linked modules or species with few partners. Since weak mutualistic selection results in weak indirect effects [16], we expect trait matching to be maximized wherever the selective pressures of interactions are more homogeneous between partners or are fewer. Thus, networks with modular structures will best allow the emergence of trait complementarity between partners despite the weak mutualistic selection. When mutualistic selection is strong, indirect effects are strong and non-interacting partners become more important in shaping coevolution [16]. Under these circumstances, indirect effects will make highly connected networks an optimal structure for the evolution of trait matching.
At the module scale, we found large poorly connected modules attained greater trait cohesion than small poorly connected ones. Yet, these results did not differ from those found in the null coevolution scenario. The similarities between scenarios suggest that comparable degrees of trait matching at a module arise regardless of coevolution. This is in agreement with the fact that the cophylogenetic signal of interacting partners is weakest at the module scale for pollination networks [32].
We detected differences in the degree of species' trait matching contingent on the role they play in the network. When coevolution is weak, we found species with both low within-module degree and AMC attained a higher degree of trait matching with their partners. Conversely, when coevolution is strong, species with both high WMD and AMC achieve greater trait matching. In other words, peripheral royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211291 species sensu Olesen et al. [30] achieve the greatest levels of trait matching when coevolution is weak, while network hubs attain the largest degree of trait matching when coevolution is strong. This change may be due to indirect effects. If indirect effects are weak when coevolution is weak, then trait evolution will be mainly driven by directly interacting species. Hence peripherals, which have few interactions and most of them are distributed in their module, will likely attain greater trait matching than other species types. When coevolution is strong and indirect effects are strong, trait evolution will be shaped by non-interacting partners. Here, super-generalists (i.e. network hubs) stand to gain since they are highly connected. Thus, hubs may not only be important in maintaining the structural coherence of networks and modules [30] but also in facilitating the matching of species traits when coevolution is strong.
When merging scales, network properties tend to outrank those of modules or species in explaining differences in trait matching between partners. Specifically, we observed the variation in trait matching at the species level to be mainly associated with the variation of structural properties at the network scale. Furthermore, the network scale is increasingly relevant as the strength of mutualistic selection increases. Though we do not directly address the question regarding the scale at which coevolution operates [1,30], our findings suggest a greater importance of network properties in constraining coevolution in mutualistic interactions.
Our work makes at least three crucial simplifying assumptions. First, we only explore how network structure influences evolutionary processes. However, the inverse relation is equally as relevant. Network structure influences evolutionary processes and, in turn, evolutionary processes shape structure. Disentangling the relative importance of each component will require models that explore both sides of the relation and their feedback.
Second, our framework does not incorporate feedbacks between ecology and evolution, which could potentially provide valuable insights into the ecological implications of the findings described here. Interestingly, de Anderazzi et al. [33] have shown that the stability of antagonistic networks increases after incorporating a feedback between trait coevolution and population dynamics. Third, our results are, to some extent, influenced by the model we used to build our networks. The niche model allowed us to build a fine-grained axis of network structure, yet it is undeniable that this model cannot fully capture the intricacies of empirical networks. In this sense, the lack of differences between the null and coevolution scenarios at the module scale could suggest that the niche model did not recreate biologically meaningful modules. The use of empirical networks or other network models may enlighten how coevolution influences trait matching at the module scale. Having stated some of our assumptions and abstractions, we argue that coevolution in mutualistic networks is neither solely driven by the strength of mutualistic selection, nor entirely constrained by the structure of the network where it is unfolding. Rather, trait matching is a joint result of selection pressures and network structures.
Trait complementarity between partners is a crucial property as it has been shown to increase the stability of mutualistic networks [18]. Human-induced perturbations are likely reshaping coevolutionary dynamics by simultaneously altering both the strength of mutualistic selection and by disassembling the web of life [34][35][36]. We expect trait complementarity between partners to decrease if perturbations reduce the strength of coevolutionary selection, as in the case of extinctions or changes in phenology. In addition, the fact that network structure can amplify the dissonance between partners further highlights the importance of conserving network structures. [37]. Yet, which architectures to preserve will depend on the strength of coevolution.