Effects of species traits, motif profiles, and environment on spatial variation in multi‐trophic antagonistic networks

Understanding drivers of antagonistic interactions across temporal and spatial scales is important for predicting community structure. In particular, studies examining spatial variation in ecological networks are critical for anticipating community responses to anthropogenic change. Most studies examining spatial interaction turnover focus on bipartite networks begging the question of whether the results are also reflected in unipartite, multi-trophic networks. To examine the spatial turnover in food web interactions, the environmental and ecological drivers of this, and the influence of interaction turnover on the preservation of individual species’ roles, we used a spatially expansive multi-trophic antagonistic ecological network data set of 129 lakes spanning over 1000 kms. We used b-diversity metrics to quantify spatial turnover in interactions and calculated the relative contributions of interaction rewiring and turnover in top, intermediate, and basal species to network turnover. We then investigated the relative and combined role of multiple ecological drivers (e.g., abundance, thermal tolerance) and environmental drivers (e.g., latitude, total phosphorus) on internal network structure. Finally, we used a motif analysis to measure the effect of spatial interaction turnover on the variation in individual species’ roles. We observed high interaction turnover across lakes, driven primarily by turnover in basal species but also the rewiring of interactions among shared species, driven, in part by underlying environmental gradients (e.g., species richness). Contrary to previous food web models applied to single sites, none of the ecological drivers we considered were effective predictors of lake-specific interactions perhaps indicating an important distinction between network model accuracy at regional and local extents. Finally, despite high spatial turnover in interactions, species’ roles were highly conserved across the study lakes demonstrating the potential of species’ roles for predicting community structure. These findings demonstrate how integrating species’ fundamental roles into trait-based approaches may improve our predictions of ecological networks at local scales.


INTRODUCTION
Ecological networks can be very complex, characterized by closely connected species and dense patterns of interactions (Montoya et al. 2006). Understanding the drivers of these interactions at multiple temporal and spatial scales is critical for predicting community structure and dynamics (e.g., Ara ujo and Luoto 2007, Morales-Castilla et al. 2015). One way of determining drivers of network variability is through measurements of network dissimilarity (e.g., Petanidou et al. 2008, Novotny 2009, Poisot et al. 2012. These analyses extend the use of b-diversity metrics (see review in Anderson et al. 2011) to investigate site-to-site variation in species' interactions (e.g., Novotny 2009, Poisot et al. 2012, CaraDonna et al. 2017. While studies of temporal turnover in interactions across seasons and years demonstrate that species interactions are ephemeral (Petanidou et al. 2008, MacLeod et al. 2016, CaraDonna et al. 2017, similar studies of spatial interaction turnover have been lacking (but see Trøjelsgaard et al. 2015, Poisot et al. 2017). Moreover, the spatial studies are focused primarily on bipartite networks involving plantpollinator (Carstensen et al. 2014, Simanonok and Burkle 2014, Trøjelsgaard et al. 2015, Emer et al. 2018, Tim oteo et al. 2018, host-parasitoid (Poisot et al. 2012, Baker et al. 2015, or consumer-resource (Novotny 2009, Lu et al. 2016, Kemp et al. 2017, Saavedra et al. 2017) systems leaving a gap, to our knowledge, in the study of spatial interaction turnover in multi-trophic, antagonistic (i.e., unipartite) networks.
The drivers of interaction turnover may not be that different from the drivers of b-diversity across large spatial or temporal gradients (e.g., Trøjelsgaard et al. 2015, CaraDonna et al. 2017 since one-way interaction turnover can occur is through a change in species composition between sites resulting in novel interactions. But the patterns may be more complex, as interaction turnover (b INT ) is the additive result of two distinct processes; differential interactions due to changes in species composition between sites (b ST ), and site-to-site differences in realized interactions between species co-occurring at sites reflecting a rewiring of interactions between sites (b RW ; Petanidou et al. 2008, Novotny 2009, Poisot et al. 2012; see Bartley et al. 2019 for further discussion). Moreover, the two sub-components of interaction turnover (i.e., b ST and b RW ) can combine to alter interactions between shared species in surprising ways. For example, in a comparison of lakes with and without invasive Micropterus dolomieu (smallmouth bass), Vander Zanden et al. (1999) found that the presence of smallmouth bass caused Salvelinus namaycush (lake trout) to shift their diet from 60% forage fish and 40% pelagic zooplankton to 20% forage fish and 80% pelagic zooplankton. Thus, a species invasion introduced new interactions between the invader and native species, but also caused native predators to shift prey resources (Vander Zanden et al. 1999). Both of these outcomes increase interaction turnover but arise from different mechanisms with different implications for forecasting the re-assembly of ecological networks under global change.
Community network structure arises through a combination of ecological drivers including random assortments of interactions proportional to species' abundances (V azquez et al. 2009, Olito and Fox 2015, CaraDonna et al. 2017, species' trait mismatch (i.e., forbidden links; Olesen et al. 2011, Ekl€ of et al. 2013, Gonz alez-Varo and Traveset 2016, and environmental constraints (Post et al. 2000). One of the more tangible consequences of global change is fluctuations in local trait distributions, potentially leading to the uncoupling of traits between species (e.g., Olesen et al. 2011, Gonz alez-Varo andTraveset 2016), a de-coupling which could be causing the rewiring of interactions across spatial gradients, potentially altering network structure and function. It has been recently demonstrated, however, that many mechanistic explanations of species' interactions are successful at predicting the summary statistics of networks (e.g., connectance), but fail to predict species-specific interactions (e.g., Allesina et al. 2008, V azquez et al. 2009, Olito and Fox 2015, Simmons et al. 2019. Studies testing mechanisms for interaction turnover rarely validate that their proposed mechanisms are correctly predicting observed interactions, and thus, the observed similarity between empirical and predictive models of interaction turnover may be arising for the wrong reasons (Olito and Fox 2015). Moreover, the implications of this interaction turnover on local food web structure and function are rarely investigated.
Quantifying the importance of individual species to community structure and function is an essential aspect of ecology, especially when predicting community response to species extinctions (e.g., Lewinsohn and Cagnolo 2012). A promising approach for measuring this preservation of species' roles across networks is based on a species' configuration within three-species subnetworks, termed motifs (e.g., Stouffer et al. 2012). These three-species motifs consist of 13 possible recurring three-node subgraphs observed across a variety of networks from biomolecules to neural cells (see Milo et al. 2002 for more discussion). A recent study by Simmons et al. (2019) used a motif approach and demonstrated that it was significantly more effective at describing network structure than traditional network metrics. Further, when Baker et al. (2015) used this approach to measure the spatial and temporal variation in individual species roles in 22 host-parasitoid communities, they observed that despite a 50% turnover in species and a 70% turnover in interactions between years, species' roles were highly conserved. This demonstrates that, in this system, the fundamental network structure was resilient to turnover in species, reassembling to preserve fundamental roles. A more recent study of spatial plant-pollinator networks in the high arctic, however, demonstrated a substantial community turnover that was related to significant variation in species' roles (Cirtwill et al. 2018). Thus, which trend holds true for the turnover in interactions across spatial or temporal gradients in multitrophic antagonistic networks remains to be tested. Consequently, a frontier in the study of ecological networks under global change is both understanding the discrepancy between the predictive accuracy of network properties (e.g., connectance, spatial turnover in interactions) and species-specific interactions (V azquez et al. 2009, Olito and Fox 2015, Simmons et al. 2019) and how individual species roles are altered by interaction turnover Baker et al. (2015). Together, they have the potential to improve the predictive capacity of our network models.
Here, we use a novel and spatially expansive freshwater fish-only food web data set (n = 129 lakes, 1094 km from south to north; Fig. 1) to quantify the spatial turnover in food web interactions, examine the environmental and ecological drivers for this turnover, and assess the implications on individual species' roles (see Fig. 2). Specifically, we examine (1) the relative contributions of interaction rewiring and turnover in top, intermediate, and basal species to overall spatial network turnover, (2) the potential environmental and ecological drivers of this turnover, including the relative and combined effect of abundance, phylogenetic relatedness, competition, thermal tolerance, and morphology on internal network structure, and (3) the implications of these changes by quantifying the fidelity of individual species' roles across networks.

Study area and data collection
Data were collected as part of the Broad-scale Monitoring Program run by the Ontario Ministry of Natural Resources and Forestry (Sandstrom et al. 2013). For this program, fish abundance and fish stomach content information was collected from 741 lakes across Ontario, Canada, between 2008 and 2012. For each lake, environmental parameters including temperature profiles (°C), lake area (ha), mean depth (m), number of degree-days above 5°C, pH, and total phosphorus (TP) were measured in conjunction with fish sampling (see Appendix S1). We used the fish stomach content information to build binary species interaction networks (where a 1 indicates the presence of an interaction and a 0, the absence) at both the individual lake and metaweb (i.e., food web of all fish-only interactions observed in stomach content of fish from the 741 lakes) level.
Studies of spatial interaction turnover are robust to sampling effects when measured against the metaweb (Poisot et al. 2012), while other studies have shown that interaction turnover between sites may be more susceptible to sampling effort (Novotny 2009). Thus, we used rarefaction techniques to minimize sampling bias (see Appendix S2). Ultimately, 129 lakes consisting of 30 different species and 1793 interactions met our criteria for inclusion. These criteria required lakes to have more than five observed interactions and a sampling effort sufficient to detect a minimum of 75% of the pairwise interactions occurring in each lake (S Obs /S Chao1 9 100; see Fig. 1 for locations and Appendix S2 for details on S Obs and S Chao1 ). In addition, we reanalyzed our data after removing the rare interactions (i.e., those interactions occurring once in a lake) to assess the influence of rare interactions on species interaction patterns (see Appendix S3 for results with rare interactions removed).

Spatial interaction turnover
To quantify spatial turnover in food web interactions, the relative contributions of interaction rewiring and turnover in top, intermediate, and basal species were calculated between all pairs of lakes (Fig. 2). Specifically, following Novotny (2009), we calculated species' interaction turnover ❖ www.esajournals.org between lake pairs using Whittaker's dissimilarity index: where a is the number of interactions shared between lakes, and b and c are the number of interactions unique to Lake 1 and Lake 2, respectively. In this way, b INT ranges from 0 (no turnover in interactions) to 1 (complete turnover in interactions). We can further differentiate b INT into its two components, turnover in species involved in the interactions (b ST ) or rewiring of interactions between shared individuals (b RW ). The contributions of b ST and b RW to b INT can be determined using similar methods laid out in Novotny (2009). In this way, both b and c, or those interactions unique to Lake 1 (b), or Lake 2 (c) can be broken down into their rewiring and species turnover components such that where b RW represents the interactions unique to Lake 1 that are due to rewiring, and b ST the interactions unique to Lake 1 that are due to species turnover (and similarly for c, those interactions unique to Lake 2). Note that this derivation of b ST and b RW differs from Poisot et al. (2012) in one critical way-our denominator (and that pre-  With this in mind, eq. (2) can be re-formulated as: And thus, interaction turnover (b INT ), can be considered the additive result of two distinct processes; (1) turnover in one or both of the species involved in the interaction (b ST ), or (2) rewiring of interactions between shared individuals (b RW ): By partitioning b INT in this way, we are able to differentiate whether turnover in species interactions across lakes is driven by turnover in species composition (b ST > b RW ), or due to rewiring of interactions between shared species (b ST < b RW ). This allows us to determine whether the cooccurrence of species is sufficient to assume an interaction (b RW = 0), or whether other environmental and species-specific factors influence turnover in species interactions between shared species (b RW > 0).
Turnover due to changes in species composition, b ST , can be further broken down (see Novotny 2009, Simanonok and Burkle 2014, Trøjelsgaard et al. 2015 where the first three b i s are turnover due to only one interacting partner being present in one of the sites (novel basal (b B ), intermediate (b I ), or top (b T ) species), and b N occurs when both interacting partners are found in only one of the sites. (1) Spatial turnover in interactions across lakes (b INT = 1 in the diagram above) is quantified using Whittaker's dissimilarity index, which is the summation of the turnover in interactions due to rewiring of interactions between species common to both sites (b RW = 0.29) and the turnover in interactions due to change in species compositions (b ST = 0.71), which can be further compartmentalized into turnover driven by changes in predator species (b T = 0), intermediate species (b I = 0.71), prey species (b B = 0), and both predator and prey species (b N = 0). Here, light gray circles with letters indicate unique species. (2) The contributions of both environmental and ecological drivers to spatial interaction turnover are calculated. (3) The implications of spatial turnover on individual species are determined by examining the fidelity of species roles as measured by their frequency in unique positions of three-species motifs, indicated by light gray circles.

Environmental drivers of spatial interaction turnover
The change in community composition over environmental gradients has been well studied (see Anderson et al. 2011 for detailed review). Few studies, however, have examined changes in species interactions over large spatial extents (Trøjelsgaard et al. 2015, Poisot et al. 2017. To determine the influence of environmental factors including latitude, area, mean depth, pH, morphoedaphic index (MEI), TP, degree-days above 5°C (DD5), and species richness on interaction turnover across lakes (see Fig. 2), we first tested for collinearity (see Appendix S4: Table S1 for the Pearson correlation coefficients). We then used a permutational multivariate analysis of variance (perMANOVA) using distance matrices to examine the relationships between each b (i.e., b INT , b RW , b ST , b B , b I , b T , and b N ) and the full set of environmental variables (i.e., b~Latitude + Surface Area + Depth + pH + MEI + TP + Richness). We conducted these analyses both with and without highly correlated environmental variables (those with a correlation coefficient greater than | 0.35|). There is no evidence to suggest an interaction between the environmental variables will improve fit (e.g., Dolson et al. 2009), and therefore, we did not include interactions among variables in our models. We used the adonis function in R package vegan package in R (version 3.5.1) and 5000 permutations (Oksanen et al. 2018; see Appendix S4 for parameter specific hypotheses).

Ecological drivers structuring interactions
Studies have shown that food web interactions can be predicted by a few trait dimensions (e.g., Ekl€ of et al. 2013). Therefore, we hypothesized that turnover in interactions across lakes could be due to ecological drivers including combinations of species abundances and trait-based explanations. We used a trait-based approach to see whether lake and species-specific drivers could be used to explain lake-level realizations of the metaweb, and, by extension, explain the large amount of interaction turnover across lakes. We tested that the following five ecological drivers structure fish interaction networks and that between-lake variations in these drivers may cause interaction turnover.
Abundance.-Neutral theory suggests that all species are ecologically equivalent and that ecological network patterns result from random interactions of individuals (e.g., V azquez et al. 2009). Therefore, we predict that more abundant species will have a higher probability of interaction than rare ones.
Phylogenetic relatedness.-It has been demonstrated that closely related species share similar traits determining trophic interactions. Further, it has been shown that the more similar species are the higher their likelihood of interaction (Cattin et al. 2004, Rezende et al. 2009). Therefore, we predict a high likelihood of interaction in closely related species.
Competition. -Species with niche overlap will compete for shared prey species and decrease the probability of interaction between a predator and a shared prey. As such, coexistence of predators with shared prey should have increased competition for prey species. Therefore, we predict that the greater the number of co-existing predators, the lower their probabilities of interacting with shared prey.
Thermal tolerance.-Lakes are divided into distinct microhabitats distinguished by speciesspecific traits, such as thermal preference (e.g., Bartomeus et al. 2016). Depending on the temperature profile of the lake, the frequency with which a predator might overlap with different prey species will change (Bartomeus et al. 2016). Therefore, we predict that greater overlap in thermal optima between predator and prey species should lead to higher probabilities of interaction.
We developed a metaweb of species interactions using all observations of interactions across the 741 lakes in Ontario for which we had stomach content data. We used data from all 741 lakes to provide the most complete representation of the possible interactions observed from the empirical data.
We then constructed lake-specific interaction probability matrices or estimates of the probability of choosing a given interaction out of all possible species interactions within a lake. First, we subset the metaweb to include only lake-specific species to identify all potential feeding interactions within a lake and then used the ecological drivers, and all possible combinations of these, to assign probabilities of observing each of these potential interactions (n = 32 combinations; Appendix S5 for more details). Probabilities for all combinations of the five base ecological matrices were calculated by ❖ www.esajournals.org multiplying the respective base model species pair probabilities together. In this way, two species with high (low) probabilities of interacting in the respective base models will have a high (low) probability of interacting in the combination of these models (see Appendix S5: Fig. S2). This allowed us to determine the relative and combined role of each driver for describing the observed species interactions within lakes. Following V azquez et al. (2009) and CaraDonna et al. (2017), we normalized the lake-specific probability matrices so that each matrix summed to one. This normalizing step was performed for all five base ecological models (e.g., abundance [A] or phylogenetic relatedness [P]) and all combinations of these base models (e.g., A9P, A9M, A9T; see Appendix S5: Fig. S2). We also created a lake-specific null model based on the metaweb where equal interaction probabilities were assigned for all species found in the lake that interact in the metaweb.
We used these lake-specific probability matrices to generate 5000 predicted interaction matrices for each lake (1 indicates the presence of an interaction, 0 the absence and the likelihood of an interaction is directly related to its probability in the lake-specific probability matrix). To do this, we fixed the number of interactions occurring at each lake to the number of interactions observed and fixed the species present at each lake to those sampled in the lake and then sampled interactions without replacement. From these lake-specific predicted interaction matrices, we calculated the specific True Skill Statistic (TSS; ranging from 1, perfect fit, to À1, inverted fit) and sensitivity value (proportions ranging from 0, poor fit, to 1, perfect fit). While we present results only for well-sampled lakes, we recognize that the absence of an interaction does not necessarily imply a true absence of that interaction in nature (see Bartomeus et al. 2016). The sensitivity values, however, allow us to avoid the problem of false absences when assessing model fit. We measure TSS and sensitivity as: where u are the interactions that were both predicted and observed to occur, v are the interactions which were predicted to occur, but observed absent, w are the interactions which were predicted to be absent, but observed to occur, and x are the interactions which were both predicted and observed to be absent. Because we are using stomach content data as evidence of interactions, we only predict interactions for those species for which we have stomach content data (i.e., piscivorous species). This consideration makes the model simulations directly comparable to the observed data for each lake. Finally, by the nature of their construction, the interaction probability matrices overestimate the prevalence of cannibalism (e.g., a species will completely overlap with its own thermal optima); thus, we removed the possibility of cannibalism in our interaction probability matrices. There were 0.13 AE 0.34 cannibal interactions per lake (mean AE standard deviation) and thus not a major contribution to network structure.

Fidelity of species roles
Finally, we determined how interaction turnover influenced species' roles (Fig. 2). Following Kashtan et al. (2004) and Stouffer et al. (2012), we used the frequency with which each species in each lake was found at each position within three-species motifs to quantify species' roles in their networks (see Baker et al. (2015) for more discussion in a bipartite network). There are 13 possible unique motifs and 30 possible unique positions within these three-species motifs (Stouffer et al. 2012). We quantified species' roles by calculating the frequency, f ij , that an individual species i, appears in each position j, within a lake normalized by that species' degree within that lake since species with a higher degree will have a larger value of f ij . In this way, each species within a network is then represented by a vector f 1 = {f i1 , . . ., f i30 }, which is henceforth referred to as the species' role.
We followed the methods outlined in Baker et al. (2015) and determined species' role conservation across lakes by performing a perMA-NOVA. Here, total dissimilarity D across all species and networks was calculated with where N is the total number of species' roles (across all species and lakes) and e ij is the distance between role i and role j. The distance ❖ www.esajournals.org metric, e ij , was calculated using the Bray-Curtis dissimilarity metric where F ij is the sum of the lesser frequencies for only those motif positions shared between species' i and species' j.
Comparisons of species' roles across networks are determined by grouping species' roles by species identity and comparing across lakes. Thus, the within species' dissimilarity for any species k is determined by where g k is the total number of species' roles across all lakes, and d ij|k is Kronecker's delta (where d ij|k = 1 if motif position k is shared between species' i and j and d ij|k = 0 otherwise).
To test for significance, we conducted two randomizations to serve as null models. The first was a true randomization where species' identity within the lakes in which they occurred was randomly shuffled. This first randomization procedure can determine whether species' roles are more conserved than random but lend little more insight into why. The second randomization assumes that species with similar generality have similar species' roles. To this end, the second randomization kept species' generality constant and shuffled species' identities within lakes with those species with the same lake-specific generality. To determine overall significance of species' roles, we compared observed D to the randomized D (D RAND and D GEN ). We repeated each randomization 5000 times. Further, we compared the reshuffled d k 's (eq. 10) to the observed d k 's to determine which species contributed more, or less, to the observed variation in species' roles. Here, we directly compared the proportion of randomized d k 's that are of equal or greater similarity than observed empirically. All analyses were performed in R 3.5.1 (R Core Team 2017) and the data and code are provided (https:// figshare.com/s/63db9179b0ec429e6b00).

Spatial interaction turnover
Overall turnover in interactions across lakes (b INT ) ranged from 0 to 1 with a mean of 0.76 (Fig. 3). We further observed that turnover in species composition (b ST , mean = 0.45) contributes more to b INT than interaction rewiring (b RW ; mean = 0.31; Fig. 3); however, there was a lot of overlap between b ST and b RW demonstrating that for some lake pairs, b RW was a more important component of b INT than b ST . Further breakdown of b ST into turnover of top, intermediate, basal, and novel components (b T , b I, b B , and b N ) indicated that turnover in basal species (mean = 0.17) is the main component of b ST followed by turnover in top species (mean = 0.14; Fig. 3).

Environmental and ecological drivers
Environmental variables explained less than 40% of the variation in overall interaction turnover (b INT ), and the variables which explained the most variation were species richness, mean depth, and latitude (R 2 = 0.12, 0.12, and 0.06, respectively; see Table 1 and Appendix S4: Fig. S1). This was primarily due to the influence Fig. 3. Comparison of each interaction turnover component where error bars represent 95% confidence intervals, boxes represent interquartile range, and middle line represents overall mean. Here, b INT represents the overall turnover in interactions between lakes, b RW the turnover in interactions between lakes due to rewiring of the interactions between co-occurring species, b ST the turnover in interactions between lakes due to species turnover, b T the turnover in interactions due to changes in top species, b I the turnover due to changes in intermediate species, b B the turnover due to changes in basal species, and b N turnover due to changes in both predator and prey species. of environmental drivers on b ST which explained almost 50% of the variation in interaction turnover due to changes in species composition. Again, this was driven primarily by depth, species richness, and latitude (R 2 = 0.17, 0.12, and 0.10, respectively; see Table 1). Rewiring (b RW ), however, was relatively unrelated to environmental drivers as they explained less than 10% of the variation, with only pH a significant predictor, (R 2 = 0.03 Table 1). The correlation analysis demonstrated that MEI and degree-days greater than 5°C were correlated with many other parameters (correlation coefficient >0.35), while depth and TP were correlated (see Appendix S4: Table S1). Even after removing MEI, degree-days, and TP from the model, the above patterns held (see Appendix S4: Table S2).
While environmental variables do predict some of the variation in interaction turnover across lakes, there is still unexplained variation. Hence, we hypothesize that the sizeable turnover in interactions across lakes (Fig. 3) could be due to ecological drivers including combinations of species abundances and trait-based explanations.
To that end, we used a trait-based approach to see whether lake and species-specific drivers could be used to explain lake-level realizations of the metaweb and, by extension, explain the large amount of interaction turnover across lakes. The five ecological drivers, abundance (A), phylogenetic relatedness (P), competition (C), thermal tolerance (T), and morphology (M), and any combination of these drivers, were poor predictors of species-specific interactions per lake (TSS < 0.23, sensitivity <0.39; Table 2). However, almost all ecological drivers performed better than the null model which had a TSS of 0.096 and a sensitivity of 0.283. The exceptions were the competition model, both alone and in combination with phylogeny and morphology. Abundance, phylogeny, morphology, and thermal tolerance combined was the best model (TSS of 0.220 and sensitivity of 0.388); however, the addition of phylogeny, morphology, and thermal tolerance only had marginal contributions to the model which incorporated only abundance (TSS of 0.196 and sensitivity of 0.365; Table 2).

Fidelity of species roles
The species-level perMANOVA for both the random null model and the generality null mode demonstrated that species' roles explained a high amount of role variability across lakes compared to the null model (observed D W ( D RAND & observed D W ( D GEN ). These results are summarized in Fig. 4, where the observed species' dissimilarity index is shown by a horizontal line, while the violin density plots are the null model results. Those species' results which are significantly different from the null model results have a horizontal line significantly lower than the violin plots and are denoted by an asterisk. When examining how individual species contribute to overall species fidelity, we observed high role Table 1. The results from a perMANOVA with distance matrices with 5000 permutations on different community interaction turnover metrics and environmental parameters. Notes: MEI, morphoedaphic index; TP, total phosphorus. Here, values represent adjusted R 2 s where bolded values indicate P < 0.05. The environmental parameters are all lake-specific values where TP stands for TP, and DD5 stands for degree-days above 5°C. The community interaction turnover metrics are defined as follows: b INT : overall interaction turnover between lakes, b RW : turnover in interactions between lakes due to rewiring of interactions between shared species, b ST : turnover in interactions between lakes due to species turnover between lakes, and sub-components of b ST , specifically b T : turnover in interactions between lakes due to novel top species, b I : turnover in interactions between lakes due to novel intermediate species , b B : turnover in interactions between lakes due to novel basal species, and b N : turnover in interactions between lakes due to both novel predator and prey species (see Appendix S3: Table S2 for the equivalent table after the removal of rare interactions).
fidelity across all species with the exception of Luxilus cornutus (common shiner), Ambloplites rupestris (rock bass), Lota lota (burbot), and Sander canadensis (sauger; Fig. 4a, c) compared to the random null model. These results differed for the generality null model where only seven species had observed species' roles significantly different from the null model (Fig. 4b, d). Overall, these analyses suggest that species' roles in freshwater fish are highly conserved across the different lakes, but that these species' roles can be predicted, for some species, by their generality.

DISCUSSION
We provide one of the first broadscale assessments of spatial interaction turnover in multitrophic antagonistic networks. We examine the relative and combined contributions of multiple ecological drivers and the environment to this turnover in a novel lake data set which spans over 1000 kms and includes over twice as many sites as the next largest spatial interaction turnover study (e.g., n = 51 in Poisot et al. (2017), n = 22 in Baker et al. (2015), and n < 20 in Carstensen et al. (2014)). Here, we find substantial interaction turnover across lakes with a larger contribution of interaction rewiring than previous bipartite food web studies (Novotny 2009, Kemp et al. 2017), but with the contribution of species turnover still as the primary driver. In addition, we find that most environmental variables were poor predictors of spatial interaction turnover and that the ecological drivers we considered were inadequate predictors of species-specific interactions per lake. Despite the substantial turnover in interactions observed across lakes, we found that species maintain high fidelity to their fundamental roles. This preservation of species' roles, especially given the large spatial extent of our study, demonstrates that despite prey switching by predators, specific network structures are maintained. In many cases, species' roles are constrained by species' generality illustrating the potential power of combining functional guilds with network models to predict local community structure.
Quantifying spatial turnover in species interactions enables us to better understand the dynamic nature of food webs across space, especially since historically, feeding interactions have been considered invariant, particularly within the same habitat type and geographic region Table 2. Results from the Monte Carlo process used to generate lake-specific interactions based on five base ecological interaction probability matrices and all subsequent combinations of those five base ecological models (see Appendix S3: Table S1 for the equivalent table after the removal of rare interactions). ‡ TSS values demonstrate overall model accuracy between the Monte Carlo process and the observed interaction matrices for each lake.
§ Standard deviations for both TSS and sensitivity are in brackets.
¶ Sensitivity assesses model fit of only interactions which are observed present.
(e.g., Morris et al. 2014). We observed a high degree of interaction turnover across lakes. The spatial interaction turnover we report (b INT = 0.76) is, on average, lower than recent studies of spatial turnover in bipartite mutualistic networks (e.g., mean b INT~0 .9 in Trøjelsgaard et al. (2015), mean~0.95 in Simanonok and Burkle (2014) and other studies of spatial turnover in bipartite antagonistic networks (e.g., mean b INT 0.85 in Novotny (2009) and mean b INT~0 .97 in Kemp et al. (2017). Consequently, current evidence suggests that overall interaction turnover may be lower in multi-trophic networks than bipartite networks, but is far from invariant as previously suggested (Morris et al. 2014). When b INT was deconstructed into b ST and b RW , b ST Fig. 4. Species-specific dissimilarity (d k ) index comparing conservation of species' roles between networks for (a & b) prey species, and (c & d) predator species where horizontal lines indicate the observed species dissimilarity index, violin plots show the results of the randomization of species roles, and (*) indicate species for which the observed species dissimilarity index was less than the randomization 5% of the time. Finally, (a & c) represent the random null model results while (b, d) represent the generality null model results. Species' names have been abbreviated with B. Stickle representing brook stickleback, C. Shiner is common shiner, J. Darter is johnny darter, N. Stickle is ninespine stickleback, R. Smelt is rainbow smelt, R. Bass is rock bass, S. Shiner is spottail shiner, W. Sucker is white sucker, N. Pike is northern pike, L. Trout is lake trout, Y. Perch is yellow perch, L. Whitefish is lake whitefish, and S. Bass is smallmouth bass. Note: Here, we classify prey species as those species for which we do not have stomach content analyses for since this is a multi-trophic system and all predators are also prey in some cases. See Appendix S3: Fig. S2 for the equivalent figure after the removal of rare interactions. contributed more to b INT than b RW (mean 0.45 and 0.31, respectively; Fig. 3) as has been shown in other studies of bipartite food webs. However, in the bipartite food webs studied by both Novotny (2009) and Kemp et al. (2017), the contributions of b RW were much lower (0.08 for both) than what we observed. This suggests that consumers in our study system demonstrate greater flexibility in prey choice than those in bipartite webs where herbivores demonstrate high resource fidelity (Kemp et al. 2017). Species turnover in plant-herbivore networks also fluctuates according to strong environmental signals (e.g., Renner and Zohner 2018). The threedimensional nature of aquatic systems, especially large bodies of water like lakes, buffers against these effects to some extent and gives species time to access suitable habitats seasonally while remaining part of the broader community. The seasonal shifts that fish may experience due to these environmental changes (e.g., changes in depth preference due to changes in thermal profile) may be at a finer scale than our sampling methods could detect. Thus, while some of the interaction turnover is interpreted as rewiring, it may be due to species turnover at a microscale.
One potential explanation for the observed turnover in interactions is the underlying environmental gradients. Interaction turnover across lakes was significantly correlated with some environmental gradients, and for most of the environmental gradients, we considered there was a significant relationship with at least one component of turnover except for surface area which was never a significant predictor. Depth was a significant predictor for both b INT and b ST (R 2 = 0.12 and 0.17, respectively; see Table 1 and Appendix S4: Fig. S1), but a poor predictor for b RW . This contrasting influence of environmental drivers on different components of interaction turnover highlights the need for caution when predicting the impact of environmental change on species interactions. For example, pH has been shown to correlate positively with many structural aspects of food webs including food web size, linkage density, and complexity (Layer et al. 2010). Here, we see that while pH is a significant predictor for turnover in interactions, it is positively correlated with turnover in interactions due to species turnover, but negatively with rewiring (see Appendix S4: Fig S1) suggesting that the positive correlation between pH and food web complexity is due primarily to new species entering the food web. For example, Kortsch et al. (2015) demonstrated that the poleward movement of generalist predators in response to a changing climate within the Barents Sea has led to significant topological changes within both the pelagic and benthic food webs. Thus, uncovering the implications of rewiring due to novel species on the stability of communities is ripe for empirical investigationparticularly as the effect of global change increases.
Beyond the influence of environmental gradients, the propensity of species to interact given trait overlap is widely recognized to be an important driver of ecological network structure (e.g., Olesen et al. 2011, Ekl€ of et al. 2013, Morales-Castilla et al. 2015. For example, freshwater fish are predominantly gape limited, consuming prey within an optimal body size in relation to their size to obtain the most energy for the lowest cost (e.g., Woodward et al. 2005). And, evidence suggests that a single trait can sufficiently reproduce ecological network structure (e.g., Petchey et al. 2008) and that only a few traits are needed to successfully predict specific species' interactions (Ekl€ of et al. 2013). Encouragingly, our comparisons between ecological traits and the null model did demonstrate that the ecological drivers improved predictions relative to a null model due to the importance of species abundance on trophic interactions. However, despite evidence suggesting the importance of the traits, we modeled in determining food web interactions (i.e., phylogenetic relatedness, e.g., Rezende et al. 2009, Stouffer et al. 2012competition, e.g., Violle et al. 2011; overlap in thermal tolerance, Bartomeus et al. 2016; and morphological matching, e.g., Woodward et al. 2005), our best model (AxMxPxT, TSS of 0.22) had lower success in predicting observed and unobserved interactions (i.e., lower TSS) than those reported elsewhere in the literature. For example, a study of the influence of predator-prey body size relationships on food web structure, Gravel et al. (2013), reported TSS values ranging from 0.13 to 0.76, depending on the food web data set used. Likewise, a study of the influence of traits on food web interactions in ground beetles, Brousseau et al. (2018), reported TSS values ranging from 0.1 to 0.65, depending on the model used. Studies attempting to reproduce ecological network structure from traits, however, typically use highly resolved networks encompassing large regional networks (e.g., Allesina et al. 2008, Ekl€ of et al. 2013). In our case, that would mean predicting the highly resolved metaweb of interactions across all of our lakes, instead of the local lake level. Traits are doubtlessly important for describing species' interactions at a local scale; however, they may need to be refined from the resolution required for the aforementioned studies. For example, the traits used in this study are typically described at a population, or lakespecific level. This is, in part, because sampling at an individual level is difficult to do in a survey of this size. However, it means that we are neglecting the contribution of intraspecific variation in these traits which has been shown to be very important in maintaining species coexistence and community structure (see review in Gonz alez-Varo and Traveset 2016).
In fact, the poor predictive ability of our traits at the local scale suggests that, at least in our study system, a few traits may not be sufficient to predict local realizations of food webs. Instead, which traits are sufficient to predict interactions at such a fine scale are likely species' and location dependent. For example, for a coldtolerant species such as Brook Trout, overlap in thermal optima may drive interactions at the edges of their range while morphology may drive interactions near the center. Moreover, the results from the environmental drivers of interaction turnover demonstrated a correlation between environmental drivers and interaction turnover suggesting that synergistic effects of both environmental and ecological drivers may be responsible for local realizations of these regional food webs. Indeed, a recent study by Joffard et al. (2019) outlined a promising technique to study such synergistic effects. Environmental and ecological perturbations are realized at a local scale, and thus, reconciling the differences in the predictive capacities of ecological network models at both local and regional scales is paramount for anticipating species-specific responses to change. For example, climate warming could shift habitat utilization by species' highly sensitive to temperature and thus substantially alter network topologies. In a recent study, Guzzo et al. (2017) demonstrated that lake trout in a pristine lake shift habitat utilization in response to yearly fluctuations in temperatures, changing prey types, and, in turn, reducing energy acquisition by this top predator.
The consequences of spatial variation in species-specific interactions on individual species are not immediately apparent. Our evidence for the preservation of fundamental species' roles across the large spatial extent of the study, despite high turnover in interactions across lakes, suggests flexibility in prey choice whereby predators select prey species which maintain stabilizing network structures (Fig. 4). This is further supported by the large contribution of basal turnover (b B ) to b ST suggesting that predators show high flexibility in prey selection. These findings are similar to those observed by Baker et al. (2015), which demonstrated the resilience of fundamental network structure despite substantial turnover in species in a host-parasitoid system and provided evidence for the reassembling of networks to preserve fundamental roles. Moreover, the preservation of species' roles suggests that the prediction of species interactions could be improved by optimizing predicted interactions based on the over and under expression of different motif profiles (e.g., Bascompte and Meli an 2005) to preserve species roles. Using a second null model based on lake-specific species' generality demonstrated that in over half the cases generality could refine predictions of species' roles. Those species which showed significant conservation of their role, however, encompassed the entire range of generality seen in our study. Examining how species' functional guild influences species' roles could provide incredibly useful in predicting local realizations of regional food webs especially if ecological networks truly have a common structural backbone (Bramon Mora et al. 2018).
We provide a comprehensive study of spatial turnover in multi-trophic antagonistic freshwater food webs. Despite high overall turnover in interactions across lakes, we demonstrate that ecological drivers are poor predictors of species-specific interactions per lake, but that species' roles are highly conserved across the study region. Our findings demonstrate that incorporating species' fundamental roles into predictive food web models may be essential to improving our predictions of ecological networks at local scales, which is instrumental to anticipating the restructuring of ecological communities as we progress in the Anthropocene.