Geographical variation of multiplex ecological networks in marine intertidal communities

ecological


INTRODUCTION
The geographic variation in environmental conditions is one of the key factors determining species distributions, contributing to community assembly and species composition (MacArthur 1972, Gaston 2000, Elith and Leathwick 2009).At local spatial scales, however, biotic interactions also play an important role in determining which species ultimately constitute an ecological community (Pearson and Dawson 2003, Thuiller et al. 2004, Wisz et al. 2013).Ecological processes such as niche partitioning or habitat selection are thus reinforced by the presence of other species and the complex landscape of competition, facilitation, and predator-prey interactions that they create (Bascompte 2009, van Dam 2009, Wisz et al. 2013).The far-reaching effects of species interactions have been demonstrated both theoretically and experimentally, but has typically only considered a small subset of "strongly interacting" species (Paine 1980, Case 1990, Menge et al. 1994, Soul e et al. 2005, Benkman et al. 2008).To what extent the network of interactions among all co-occurring species varies across geographical environmental gradients has scarcely been assessed.
Changes in species diversity and composition along environmental gradients are widely recognized and well established (Gaston 2000, Hawkins and Felizola Diniz-Filho 2004, Hillebrand 2004, Fine 2015).Although it has been demonstrated that these changes in species composition across space, through mechanisms such as environmental filtering, strongly influence network structure via the addition and disappearance of species (Tylianakis andMorris 2017, Pellissier et al. 2018), we are just starting to get a picture of how bipartite interaction networks (e.g., plant-pollinator or host-parasitoid networks; Kissling and Schleuning 2015, Morris et al. 2015, Trøjelsgaard et al. 2015, Pellissier et al. 2018, Galiana et al. 2019) and food webs (Romanuk et al. 2006, Baiser et al. 2012, Cirtwill et al. 2015, Wood et al. 2015, Gilarranz et al. 2016, Kortsch et al. 2019) change across geographical gradients of heterogeneous environmental conditions.Ecological communities are however comprised by a multiplicity of interaction types (Pocock et al. 2012, Garc ıa-Callejas et al. 2018a), which joint consideration can deeply influence the way we understand ecological communities (K efi et al. 2012, 2016, Mougi and Kondoh 2012, Lurgi et al. 2016, Garc ıa-Callejas et al. 2018b).Hence, incorporating the environmental sensitivity of different types of species interactions, trophic and non-trophic, is fundamental to fully understand the geographical variation of ecological communities.
Environmental conditions can affect ecological interactions differently depending on the type of interaction.For instance, while trophic (i.e., predator-prey) and non-trophic negative interactions (e.g., competition) have been shown to be stronger (i.e., higher predation pressure and stronger competition) in warmer and more constant environments (LaManna et al. 2017, Roslin et al. 2017), facilitative interactions are expected to be more prevalent in harsh environments (Choler et al. 2001, Maestre et al. 2005, Gallien et al. 2018).Moreover, several non-trophic interactions increase the establishment of consumers that otherwise would not be able to persist in local communities.These consumers can in turn affect the facilitator species by consumption, thus creating complex dependencies between different interaction types (K efi et al. 2015, 2016).Therefore, the structure of networks of ecological interactions comprised of different interaction types is expected to change differently across environmental gradients.Several environmental properties have been related to different aspects of network structure.Environmental constancy, for example, has been shown to play an important role in structuring local host-parasitoid communities across the globe (Galiana et al. 2019), and has been related to the specialization of mutualistic networks (Schleuning et al. 2012, Trøjelsgaard and Olesen 2013, Dalsgaard et al. 2017).Yet, conclusions on the variation of the structure of ecological networks across a given environmental gradient are rarely drawn due to the idiosyncrasy of the results (Moles and Ollerton 2016).Many aspects of the analyses influence the results observed, such as the type of network properties analysed (Dalsgaard et al. 2017, Galiana et al. 2019), the spatial scale under consideration (Wood et al. 2015, Galiana et al. 2019), or the sampling effort performed (Morris et al. 2014, Fr€ und et al. 2016, Vizentin-Bugoni et al. 2016).
One fundamental aspect to consider when comparing the structure of complex ecological networks, whether across environmental gradients or in other contexts, is species richness.Several aspects of network structure are influenced by changes in species richness and connectance (Bengtsson 1994, Vermaat et al. 2009, Riede et al. 2010, Dunne et al. 2013, Poisot and Gravel 2014).For instance, the average number of links per species increases as communities become larger (Riede et al. 2010, Dunne et al. 2013).Thus, studies addressing the variation in network properties along any environmental gradient in which a species diversity gradient exists (e.g., latitudinal gradient in species richness) should account for the effect of changes in species richness on other network properties to properly unravel the effects of environmental conditions on network structure (Cirtwill et al. 2015, Galiana et al. 2019).
Here, we analyze the geographical variation of multiplex species interaction networks (i.e., incorporating several interaction types) in rocky shore marine intertidal communities across the central coast of Chile.These interactions, and the multiplex networks they create, were identified in detail in previous studies (K efi et al. 2015, 2016).In this system, geographic variation in diversity and community composition has been related to different spatially structured environmental factors, including recruitment, nutrient loadings, and thermal regimes (Navarrete et al. 2005, Wieters et al. 2009, Broitman et al. 2011).In particular, previous studies have shown that spatial variation in coastal upwelling (i.e., an oceanographic phenomenon through which cooler, and usually nutrientrich, water is transported from deep areas toward the surface of the ocean) along this region represent an important driver of observed geographical differences among rocky shore intertidal communities (Broitman et al. 2001, Navarrete et al. 2005, Valdivia et al. 2013, 2015).Coastal upwelling events are linked to increased productivity along specific regions of the coast (Tapia et al. 2014), influencing nutrient and resource availability for coastal benthic communities (Nielsen and Navarrete 2004, Navarrete et al. 2005, Wieters 2005).Productivity has been linked to different aspects of food web structure such as food chain length, the average number of prey attacked by top predator species in terrestrial semiarid systems (Arim et al. 2003, Arim andJaksic 2005), or to the relative dominance of particular species or functional groups in marine rocky shore communities (Nielsen andNavarrete 2004, Wieters 2005).More productive environments usually promote complexity in the structure of model networks of ecological interactions (Neutel et al. 2007).
Building on this knowledge, we hypothesize that environmental factors associated with upwelling events, which translate into the sustained presence of cold water nearshore, are positively correlated with the complexity of ecological networks in marine rocky shore intertidal ecosystems.To test this hypothesis, we analyze the relationship between different environmental factors, related to water temperature and oceanographic regimes, and the structure of networks composed of three different types of biotic interactions: trophic (i.e., feeding), nontrophic negative (e.g., competitive), and non-trophic positive (e.g., facilitative), that occur among the species in these communities.This study constitutes a first attempt at gaining a better understanding of the effects of geographical variation in environmental conditions on species composition and the structure of multiplex ecological networks in marine intertidal communities.This knowledge is fundamental to be able to predict the potential consequences of environmental change on broad-scale geographical patterns of variability in ecological communities.

Study location and spatial structure of local communities
To assess the geographic variability in the network of ecological interactions of the rocky shore intertidal communities along the central coast of Chile, we used data on community composition from 17 sites spread along the coast, from 29.5°S to 36.07°S (i.e., ~970 km of the shoreline), and surveyed multiple times over multiple years as part of a previous study (Broitman et al. 2011; Fig. 1 and Appendix S1: Table S1).Presence and abundance of all sessile and mobile benthic organisms were recorded at each sampling location (see Appendix S1: Section S1 and (Broitman et al. 2001(Broitman et al. , 2011) ) for additional details on sampling protocols and procedures).The sampling design was sufficient to capture robust local-scale patterns of species richness and abundance (Broitman et al. 2011).
Sampled communities comprise gastropods, bivalves, barnacles, asteroids, and decapod crustaceans, along with chlorophycean, rhodophycean, and phaeophycean macroalgae (Oliva and Castilla 1992, Hoffmann and Santelices 1997, Espoz et al. 2004, Broitman et al. 2011).The taxonomy of sampled species was resolved to the lowest possible level, with a few species such as turfforming algae recorded as genera due to the inability to classify them to the species level in the field.Hence, species in our communities represent a mixture of taxonomic resolutions, and we define species richness of communities as the total number of distinct taxonomic entities present in them.We do not expect this to influence our results in terms of community or network structure (Williams andMartinez 2000, Underwood andChapman 2006).A full list of the species considered in this study is presented in Appendix S1: Table S2.

Network construction
Multiplex ecological networks of species interactions between all the species in the data set were obtained from (K efi et al. 2015, 2016).These studies describe the whole set of ecological interactions known to occur between the species found in rocky intertidal communities of the central coast of Chile.The data set includes trophic, non-trophic negative (e.g., competition for space, predator interference) and non-trophic positive (e.g., recruitment facilitation, habitat provisioning) ecological interactions between pairs of species.From these interactions, three layers of the multiplex ecological network were constructed at the regional scale; one for each of the three different aforementioned interaction types.This regional network constitutes a potential web of interactions between species (also sometimes referred to as "metaweb").The resulting potential trophic (i.e., food web), non-trophic negative and non-trophic positive layers of the multiplex ecological network, comprising 106, 76, and 69 species, and 1,362, 3,089, and 172 links, respectively, are presented in Data S1.
Local-scale information on the species' presence/absence at each study site (as described above) was combined with the metaweb to assemble local, realized webs of interactions.All samples at each given location were pooled across years to create a single network per site.Three layers of the multiplex ecological network were defined for each of the 17 sites studied in this way; one for each interaction type.We assumed that, wherever a pair of species sharing a link in the metaweb co-occurred at any given site, they would interact.Links were thus added to local/realized webs based on species co-occurrences.Plots of the networks for all sites and interactions types can be found in Appendices S2-S4 and can be visualized interactively on an online tool for network visualization (available online).7

Probabilistic interaction networks
In addition to determining the structure of local multiplex interaction networks, we used species co-occurrences to define local interaction probabilities between species as the product of the probability of occurrence of each species involved in a given interaction at the site level (Appendix S1: Section S1).We are aware of the many limitations of the co-occurrence approach to inferring ecological interactions (Freilich et al. 2018) and their strength (Barner et al. 2018), which in some cases can even suggest opposite trends in the sign and magnitude of specific interactions (Thurman et al. 2019).In spite of these caveats, we believe that given the data we have, this is the best way in which we can derive an estimate of interaction probability.Even if it is not perfect, results based on this approach can still inform our expectations on the variability of weighted interaction networks, and results from these exercises may be used for future comparisons.This approach allowed us to build, in addition to the qualitative (unweighted) versions of the networks as highlighted above, probabilistic networks of ecological interactions with the same topological structure as their qualitative counterparts.We then statistically compare the resulting network properties of each network type across sites to examine the influence of considering more quantitative information when building ecological networks (see Statistical analyses).

Network structure
Once networks of ecological interactions were constructed at the local scale for each site, their structure was analysed.For the three interaction networks, we recorded the number of species in the network (S), the number of links (L), connectance (i.e., the fraction of realized out of all possible links, C), the average number of links per species (L/S), the average number of incoming and outgoing links per species (indegree and outdegree, respectively), and network modularity (Q; see Appendix S1: Section S1 for details on modularity calculation).The average number of incoming and outgoing links per species can be interpreted in different ways depending on the interaction type considered.In the trophic layer of our networks, those measures quantify the average number of resources (generality; Gen) and the average number of consumers (vulnerability; Vul) of the species in the network, respectively.In the nontrophic negative and non-trophic positive interactions layers of the networks, these quantities can be related to competition pressure and facilitation provided from vs. delivered upon other species, respectively.These measures provide information on how interactions are distributed across species.We additionally calculated the  S1, where more information about each site, including their full name, can be found.The structure of the trophic, non-trophic positive, and non-trophic negative layers of the network of ecological interactions is shown for the northernmost (i.e., TEM), a middle-range (MONT), and the southernmost (i.e., BUP) sites considered in this study.Species (i.e., nodes in the network) are colored according to the module they belong to extracted from the modularity analysis.
variability in generality and vulnerability values across species in the network, normalized by L/S (SD Gen and SD Vul, respectively).These were quantified as the standard deviation across individual species values, following (Williams and Martinez 2000).
In addition to the features quantified across the three network layers considered, we analysed measures that are specific to particular interaction types.In food web analyses, the mean food chain length (MFCL) is commonly used to assess the vertical diversity of communities (Vermaat et al. 2009), with more complex food webs exhibiting longer chains of contiguous links from basal resources to top consumers.MFCL is quantified by measuring the length (i.e., counting the number of links) of all the possible trophic chains from each basal resource to each top consumer and averaging across.Analysis of the vertical diversity of the trophic layer of the networks was complemented with the fraction of basal (B), intermediate (I), and top consumer species (T).These measures indicate how species are distributed across trophic levels.
For non-trophic negative and positive interaction networks, we additionally quantified the fraction of mutually competitive vs.amensalistic interactions, and the fraction of mutually facilitative vs. commensalistic interactions, respectively.Mutually competitive interactions occur when a negative interaction between two species goes in both directions, i.e., species A negatively affects species B and vice versa; whereas amensalistic interactions are those in which the negative effect only goes from one species to the other, i.e., either species A negatively affects species B or the inverse is true, but not both.Similarly, in positive interactions, a mutually facilitative interaction is one (or two independent interactions that are both beneficial) that benefits both partners, whereas a commensalistic interaction only benefits one of the species embarking on the interaction while the other species is not affected.These measures assess the degree of asymmetry of interactions in these networks.
Probabilistic interaction networks were analysed using quantitative descriptors of some of the properties introduced above.See Appendix S1: Section S1 for details on quantitative descriptors.

Environmental variables
The investigation of the effects of environmental factors on the variability of ecological communities relies on a careful selection of the appropriate aspects of the environment thought to drive the observed patterns.When hypotheses on the relationships between environment and community structure are clear, and environmental factors are well defined, a few predictive variables can provide relevant information on the environment and its effects on communities (Kortsch et al. 2019).Based on the hypothesis that, in the rocky shore intertidal of the central coast of Chile, some of the geographical variation in community structure might be driven by differences in the supply of nutrients and particles (e.g., larvae, phytoplankton) usually associated to the presence of cold water advected to the surface by coastal upwelling, we focused on a set of variables that are able to capture local upwelling dynamics.
We centered our analysis of environmental features around sea surface temperature (SST) variability over different timescales, and direct thermal indices of upwelling intensity.Variability in SST over different timescales reflects the influence of different processes.For instance, upwelling variability leaves characteristic patterns of variation in SST over synoptic timescales (Tapia et al. 2009(Tapia et al. , 2014)).The annual variability in SST is related to the predictable seasonal changes in solar radiation, and intra-seasonal SST fluctuations can be attributed to Kelvin waves (Tapia et al. 2014).Table 1 shows all the variables considered for our analysis.Measurement of these variables was based on satellite image analyses.See Appendix S1: Section S1 for full details on data sources and how the variables were measured and calculated.

Statistical analyses
Since all the predictor environmental variables described above are related to SST, it is expected that they show strong correlations with one another and in particular with the long-term average SST (LT SST; Appendix S1: Fig. S1).To reduce multicollinearity among our environmental predictors, we performed linear regressions of each of the remaining five variables against LT SST and used the standardized (i.e., scaled to mean = 0 and variance = 1) residuals from those regressions as predictor variables in subsequent analysis; except for LT SST, which remained as originally extracted.This resulted in smaller covariation among variables.LT SST has been shown before to be a robust indicator of the role of environmental variability on ecological processes across coastal upwelling ecosystems (Broitman et al. 2001, Menge et al. 2003, Navarrete et al. 2005, Wieters et al. 2009, Valdivia et al. 2013).This generalization seems to hold well along our study region (Lara et al. 2019).
To understand the effects of geographically varying abiotic conditions on network structure, while at the same time disentangling which of these effects are manifested indirectly through changes in species richness (S) from those occurring directly due to changes in these conditions, we used piecewise structural equation models (SEMs).SEMs were built by defining network properties as linear combinations of the environmental variables described above plus species richness, while species richness itself was defined as a function of the environmental variables alone.Because network properties should be allowed to freely covary across the geographical gradient we decided to adopt a similar approach to that taken in (Byrnes et al. 2011) and fit single individual SEMs for each of the network properties analysed, except for S, which featured in all SEMs.Continuous variables were fit using linear regressions whereas count data (i.e., species richness and number of links) were fitted using generalized linear models with quasipoisson and negative binomial responses to account for under-and over-dispersion in the data, respectively.Assumptions of the linear models were quantitatively validated using the gvlma package in R (R Core Team 2017) and, occasionally, the Breusch-Pagan test for heteroscedasticity.All linear regressions performed complied with the assumptions.Goodness of fit of generalized linear models was quantified using the significance of the similarity of the models' deviance to a chi square distribution.Further details on the specification of SEMs can be found in Appendix S1: Section S1.
To assess the relationship between the variability in qualitative networks vs. the variability in the probabilistic ones, we used Procrustes analysis, testing for its nonrandomness using the protest statistical test (see Appendix S1: Section S1 for full details).
The data used in this study, and the source code used to analyze them, including all statistical analyses, can be found in Data S1.All statistical analyses were performed in R (R Core Team 2017).SEMs were fitted using the piecewiseSEM package for R (Lefcheck 2016).

Trophic interactions
The structure of the trophic layer of the multiplex network of ecological interactions in rocky shore communities varies considerably along the central coast of Chile (Fig. 1).Among the environmental variables analysed, the best predictors of food web structure were: (1) the long-term average of the sea surface temperature (LT SST), (2) the fraction of days of the upwelling season during which a community is exposed to temperatures below 14°C (Fr Days), (3) environmental unpredictability, measured as the fraction of variability in SST across the year that cannot be explained by the seasonal cycle (Fr Annual); (4) and the strength of the upwelling regime along the coast (Upwelling; Fig. 2a, Appendix S1: Table S3, and Appendix S1: Fig. S5).We note that the overall R 2 values for the models included in the SEMs are unusually high, but careful inspection of each component of the SEMs showed no errors in their calculations.We attribute this to the strong relationship between number of species (S) and number of links (L) and other food web attributes (Appendix S1: Figs.S5-S7).Removing S from the models as a predictor variable causes a drastic drop in all fits, but doing so also removes one of the predictor variables we want to examine.
Fr Days was positively related to the number of species (S) and the connectivity of species within the network (L/ S, Vul, and SD Vul).It was also negatively related to the mean and variability of the diet breadth of consumer species (Gen and SD Gen, respectively), the fraction of basal resources (B) and modularity (Q; Fig. 2a and Appendix S1: Table S3).LT SST seems to only affect the average number of interactions per species, modularity and the fraction of intermediate species; whereas environmental unpredictability (i.e., Fr Annual) negatively affects the variability in the number of resources across species (SD Gen), modularity, and the fraction of basal species, while being positively related to different aspects of connectivity such as L/S, Vul, and SD Vul (Fig. 2a and Appendix S1: Table S3).Finally, cold water upwelling intensity (Upwelling) seems to negatively affect network connectivity (L and L/S), variability in prey vulnerability, and the fraction of intermediate species; while favoring basal species and the variability in diet breadth across the network (Fig. 2a and Appendix S1: Table S3).
Interestingly, we found that even though some of the changes observed in food web structure across the geographical extent studied are manifested via variation in species richness (S), such as in the average number of links per species (L/S), consumers' diet breadth (Gen) and its variability (SD Gen), or the fraction of basal and intermediate species (B and I, respectively); most of these changes are independent from changes in S (Fig. 2a and Appendix S1: Table S3).These results thus suggest that the environmental variables considered here are good predictors of network structure along the geographical gradient studied and affect food web properties independently from their effect on species richness.

Non-trophic positive interactions
The structure of the non-trophic positive interactions layer of the network in the marine rocky intertidal of  S4, and Appendix S1: Fig. S6).Only the fraction of days during upwelling seasons with SSTs below 14°C (Fr Days) impacts noticeably a few aspects of network structure (Fig. 2b and Appendix S1: Table S4), such as the number of species (S), and modularity (Q; Fig. 2b).The geographical variability observed in the structure of these networks in terms of number of links (L) and the average number of species that receive positive interactions (Gen) is mediated by changes in species richness (S) (Fig. 2b and Appendix S1: Table S4).

Non-trophic negative interactions
Among the three network layers analysed, the structure of the non-trophic negative interactions layer of the network seems to be the most affected by abiotic factors (Figs. 1 and 2c, Appendix S1: Table S5, and Appendix S1: Fig. S7).The fraction of days during the upwelling seasons with temperatures below 14°C (Fr Days) was found to be positively correlated with species richness (S) and with the variability in the number of negative interactions a given species is subjected to, and exerts over, other species (SD Vul and SD Gen, respectively).Temperature (LT SST), on the other hand, seemed to have a stronger influence on network modularity and connectivity (L, L/S, Gen, and Vul; Fig. 2c, Appendix S1: Table S5).Communities exposed to lower temperatures thus comprise larger numbers of competitive interactions and to display a more modular organization.
A deeper analysis of these relationships reveals that while the effect of LT SST on network structure is very much independent from the number of species in the network (S), the effect of Fr Days on network structure is largely mediated by S (Fig. 2c and Appendix S1: Table S5).In particular, observed effects of Fr Days on the total number of links (L), and on the average number of interactions per species (L/S, Gen, and Vul) occur mainly through S (compare the significance, i.e., arrow width in Fig. 2c, between direct and indirect effects; Appendix S1: Table S5).This result suggests that while species richness is mainly related to the fraction of days during which communities are exposed to nutrient-rich waters, the organization of the interactions among these species is mostly related to the long-term average SST.
Environmental unpredictability (both in terms of Fr Annual and Clim) was found to be correlated with only a few network properties (Fig. 2c), and not as significantly by comparison.This suggests that the main influences of the environment on network properties across the geographical gradient are exerted by Fr Days and LT SST.
Are probabilistic interaction networks structured in a similar way to qualitative networks?
Our results showed that quantitative descriptors of food web properties, which incorporate information from species co-occurrence data as a proxy for interaction probability, vary in synchrony with their qualitative counterparts (i.e., those considering only the binary network data) across the study region (Appendix S1: Section S2).This was true across the three interaction types studied here.These results however, should be taken with caution, since, as mentioned above, there are many potential problems with building probabilistic networks.Our results nonetheless provide a picture of what can be expected from probabilistic ecological networks built in this way.

DISCUSSION
Our results show that the structure of the multiplex network of biotic interactions (i.e., a network including three types of interactions between the species of a community) in the marine rocky-shore intertidal of central Chile is systematically and differentially affected by environmental variation depending on the interaction type considered.The structural organization of non-trophic negative interactions is more sensitive to environmental variation than that of feeding relationships, with the structure of non-trophic positive interactions being the least labile to environmental fluctuations.We also show that environmental effects on these multiplex ecological network are partly mediated through changes in species richness and partly via direct influences on the layout of how species interact, probably associated to changes in "environmental predictability," and to bottom-up nutrient availability.
intertidal communities are affected by environmental factors.The figure shows the summarized results of saturated individual SEMs built for each network property (see Methods).For clarity, only those predictors and relationships identified as significant by SEMs (P < 0.05) are presented.Complete SEM diagrams including the full set of variables used and non-significant paths are presented in Appendix S1: Figs.S8-S10.Standardized coefficients are shown overlaid on the corresponding arrows.Red and black arrows indicate negative and positive effects respectively.Adjusted R 2 values of the fit for each model are presented inside the box of the corresponding response variable.Connectance (C) was removed from the diagrams to make them clearer but its position in the network is comparable to links per species (L/S) as a related measure of connectivity.See Appendix S1: Tables S6-S8 for the standardized coefficients of each predictor on the response variables.Q, network modularity; B, fraction of basal species; I, fraction of intermediate species.

(FIG. 2. Continued)
Communities subject to colder temperatures for extended periods of time, i.e., those where nutrients are steadily supplied from upwelled waters, tend to be more complex, both in terms of species and trophic interactions.Consumers tend to be more diet specialists with increased nutrient availability, with prey resources experiencing greater predator pressure.This is in turn associated with stronger homogeneity of trophic interactions across resource species in more persistently cold environments.Communities experiencing stronger upwelling events, on the other hand, tend to harbor smaller number of ecological interactions.In contrast, the structure of the non-trophic positive interaction network layer remains virtually unaffected by these environmental variables, being only weakly influenced by the fraction of days during upwelling seasons during which communities are exposed to temperatures below 14°C (Fr Days).Last, the structure of non-trophic negative interaction networks is more strongly related to long-term sea surface temperature (LT SST) than to the dynamical behavior of temperature regimes."Colder" communities are comprised of a larger number of negative interactions relative to their "warmer" counterparts; with the number of negatively interacting partners being more variable among species in warmer environments.
The strength and prevalence of competition, the dominant non-trophic negative interaction in our network, has been associated with environmental conditions before.Intraspecific competition between forest trees has been shown to be stronger in tropical forests relative to their temperate counterparts (LaManna et al. 2017); pointing to the role of competition in maintaining biodiversity and the potential for warmer and more predictable environments to harbor higher competition levels.This phenomenon does not seem to be exclusive to intraspecific interactions, but also shared with predator-prey interactions such as arthropods attacking herbivore prey (i.e., caterpillars), in which predation pressure increases toward the tropics (Roslin et al. 2017), where climate is warmer and more constant (see also Vermeij 1994, Freestone et al. 2011).Even though our admittedly crude measure of "environmental unpredictability" did not show up strongly in our analysis of the geographical variation of non-trophic negative interactions, we found that it strongly influences the structure of food webs in the Chilean marine intertidal community by increasing the average number of links per species and, therefore, network connectivity.
The percentage of variance in SST not explained by the seasonal cycle, termed environmental unpredictability here, has been alternatively interpreted as a measure of attenuation of the annual temperature cycle.Since upwelling of cold water from deep layers of the ocean is often strongly intensified (due to alteration of winds) during spring-summer (at least south of 32°latitude), this can be reflected in a more attenuated annual cycle in SST, which could be as predictable as when the annual cycle explains most of the SST variability.A thermal regime less dominated by the annual cycle can be considered more temporally "heterogeneous" in the scales of dynamics that influence it.Such increased temporal heterogeneity in the environment may represent a mechanism to increased species richness and/or trait tradeoffs (Blanchette et al. 2009), which can in turn influence the structure of the network of ecological interactions.
Previous studies investigating the geographical variation of the intertidal communities analysed here revealed strong correlations between environmental factors and mesoscale oceanographic processes, such as sea surface temperature and the intensity or frequency of upwelling events, and community structure.In particular, a negative latitudinal gradient was observed in the number of algal species across the central coast of Chile (Broitman et al. 2011), and, more generally, of algal and mollusk species diversity across the entire southeast Pacific coast (Santelices andMarquet 1998, Valdovinos et al. 2003).Furthermore, sea surface temperature and upwelling events are known to influence algal cover and the abundance of both sessile (including algae) and mobile organisms across this region (Broitman et al. 2001, Nielsen and Navarrete 2004, Tapia et al. 2014).Our results are in line with these previous studies and extend those analyses to different aspects of the structure of multiplex networks of ecological interactions among species.
Our findings additionally reveal that changes in species composition due to environmental factors affect these communities differently depending on the interaction type considered.Sea surface temperature remained strongly correlated with the structure of the non-trophic negative, but not the trophic or the non-trophic positive, layers of the multiplex networks of ecological interactions.On the other hand, even though our measure of upwelling intensity per se (UI) did not feature prominently in our analysis (except for a weak influence on food web structure), as opposed to previous analyses of community structure (Broitman et al. 2001), an intimately related measure, namely the fraction of days during upwelling seasons during which communities are exposed to temperatures below 14°C (i.e., water that is not nutrient-depleted (Tapia et al. 2014)) was strongly correlated with several aspects of multiplex networks of ecological interactions.This suggests that it is the length of periods of time of exposure to colder temperatures, or better, the fraction of time exposed to nutrient-rich (or at least non-depleted) waters (Tapia et al. 2014), that strongly influences community structure through bottom-up nutrient loadings that can be captured by benthic macroalgae (see Nielsen andNavarrete 2004, Wieters 2005).
Coastal upwelling events are linked to increased productivity along specific sections of the coast.These events thus affect the nutrient and resource availability of these communities, which are exposed to these events for variable periods of time across the year (Wieters 2005, Wieters et al. 2009, Tapia et al. 2014, Valdivia et al. 2015).Effects of this exposure to resource availability have the potential to propagate up the trophic and non-trophic networks of interactions to varying degrees (Nielsen and Navarrete 2004, Aldana et al. 2017, Fuentes et al. 2017).Several studies on Chilean and other upwelling rocky shore ecosystems have shown that such upwelling-driven bottom-up effects can affect the relative importance of facilitation and the intensity of top-down regulation over meso (among sites in our study) and regional scales (Menge 1992, Menge et al. 2003, 2004, Navarrete et al. 2005, Wieters 2005), and are therefore expected to affect network structure.
Our results extend previous findings by revealing that rocky shore intertidal communities less exposed to nutrient-depleted waters (temperatures above 14°C) harbor more connected food webs, in which the fraction of intermediate consumers is larger than in warmer, possibly less productive environments, increasing thus the vertical diversity of these networks.Surprisingly, we did not find strong correlations between coastal upwelling, and hence potential sources of productivity, and any aspect of the structure of competitive networks.It is possible that, considering the overall high levels of nutrients in the Humboldt ecosystem, even in weak upwelling areas, population growth rates are always sufficient to engage in competitive interactions for space.It also suggests a potential reason for a general lack of evidence for productivity-competitive networks relationship in the literature, except when looking at the indirect kinds such as apparent competition (Tylianakis and Morris 2017).
The dichotomy between the quantitative and qualitative (i.e., binary) structure of networks of ecological interactions has highlighted the necessity of considering both aspects of network structure to fully understand community organization (Bersier et al. 2002, Gilarranz et al. 2011).We do not have estimates of strength for all these species interactions, but using interaction probability as a proxy to quantify interactions, we have shown that environmental factors do not differentially affect the structure of qualitative vs. probabilistic interaction networks in our system.This congruent change in network structure across the environmental gradient studied sheds confidence on the potential interactions comprising our networks as good representations of realized ones.Admittedly, our networks of interaction probabilities are not able to capture the strengths of the interactions analysed.Nonetheless, interaction probability is a common approach to quantify interactions in ecological networks, and it is sometimes the best that can be obtained given the difficulties in collecting empirical data for the strength of ecological interactions (Poisot et al. 2016).Further studies are needed to obtain estimates of interaction strengths in order to assess whether our findings on probabilistic networks can be extended to networks of quantified interactions strengths.
Similarly, an important aspect of the present study that might potentially affect patterns of network variation across gradients, is the way species interactions were characterized.Here, species interactions were defined based on expert knowledge of the system along with considerable experimental evidence, thus rendering the multiplex network obtained a collection of potential interactions.Additionally, local network structure was inferred based on species co-occurrence at each site, which has been demonstrated to bias the inference of ecological interactions, since the realization of these might be dependent on other processes such as densitydependence or behavioral shifts (Novak 2013).It should be noted, however, that in most cases, predators' functional responses produce rapid switches in the abundance of consumed prey, but they rarely drop a prey species altogether, when the species is still present (as shown by Novak 2013; Fig. 4).Therefore, while our topological networks would not capture such dynamical and ecologically important functional responses and switches in attack rates, the approach would correctly identify the links (strong and weak) as present.This has the potential to affect network structure given that interactions considered might not actually be realized, at least at local scales.Future research aimed at improving the network of ecological interactions by collecting empirical observations of these interactions should prove useful in supporting the patterns currently observed.Patterns of local species abundance could be used to reconstruct realized interaction networks at local scales.Such measures of abundance must be expressed in a common currency (e.g., biomass) for them to be useful.This represents a major challenge when communities include sessile and mobile species, individual and colonial organisms, and/or a mixture of animal and algae species.Further theoretical developments are also needed to define how species abundances can be used to infer interactions and their strengths (Berlow et al. 2004).
Here, we focus on environmental variables related to thermal variability.These have been found to be among the main drivers of community structure in rocky shore ecosystems along the coast of Chile (Nielsen and Navarrete 2004, Wieters 2005, Wieters et al. 2009, Valdivia et al. 2013).The potential effect of wave exposure has been minimized in our studies by selecting similarly wave exposed platforms.However, other factors, which covary with SST, such as levels of nutrients, dissolved oxygen, pH, salinity, or alkalinity have the potential to affect the assembly of these communities.Unfortunately, observations for these variables across our study sites and the years of study are not available.This information could allow us to advance toward a more mechanistic model for community assembly of rocky shores.To a greater or lesser extent, all these variables are related to coastal upwelling at some spatial and temporal scales.Further work should thus focus on devising experiments and large-scale observations to find ways of teasing apart these potential driving factors from the more phenomenological driver we call "upwelling" and its associated thermal signal.This will ultimately result in a better understanding of the effects of geographical variability of environmental conditions on the structure of biodiversity.

CONCLUSION
We have shown that environmental conditions, by being strongly related to several aspects of network structure, have the potential to drive the assembly of multiplex networks of ecological interactions.Some of the changes observed in network structure across the environmental gradient are mediated by species richness, but many environmental effects are direct and strong on the networks of interactions themselves.Our work highlights the importance of considering multiple interaction types when investigating the geographical variation of ecological communities, as these respond differently to natural environmental variation, especially so to factors producing variation in productivity.Understanding these patterns of variation is not enough, however.We need network theoretical studies that are able to help us elucidate the mechanisms through which environmental factors generate the whole-community patterns observed.Future research efforts should focus on this.After all, "Understanding patterns in terms of the processes that produce them is the essence of science" (Levin 1992).Moreover, this understanding is compulsory if we are to predict the responses of these complex ecosystems to climate change.
FIG. 1. Spatial location of sampling sites and mean sea surface temperature (SST) along the central coast of Chile.Sea surface temperature (SST; °C), was obtained as the average value across science quality, global, nightly images captured by Pathfinder V5.2 at a 4-km spatial resolution (see Methods) from 1 August 1997 to 30 March 2004.Labels of the sampling locations correspond to the site codes presented in Appendix S1: TableS1, where more information about each site, including their full name, can be found.The structure of the trophic, non-trophic positive, and non-trophic negative layers of the network of ecological interactions is shown for the northernmost (i.e., TEM), a middle-range (MONT), and the southernmost (i.e., BUP) sites considered in this study.Species (i.e., nodes in the network) are colored according to the module they belong to extracted from the modularity analysis.
FIG. 2. Pathways through which geographical environmental variability affects the structure of the multiplex network of ecological interactions.Structural Equation Models (SEMs) diagrams showing the significant pathways through which properties of the (a) trophic, (b) non-trophic positive, and (c) non-trophic negative layers of the multiplex interaction network rocky shore marine

TABLE 1 .
Environmental variables used to quantify environmental geographical variability in this study.