General model
Whole food webs can be decomposed into collections of \(p\)-species interactions called motifs35 that provide a mesoscale characterization of the structural properties of ecological networks28,36–39. In a \(n\)-species food web (\(n\ge p\)), the collection of \(p\)-species motifs (\(p\le n\)) in which species \(i\) is involved in (\({M}_{i}=\{{m}_{i,1},{m}_{i,2},...,{m}_{i,x}\}\)) forms its motif census (\({M}_{i}\)) 14,29. The motif census provides an overview of all the interactions and connected species likely to affect a species’ dynamics, and the propagation of disturbances through their interactions. Here, we focus exclusively on the most abundant 3-species motifs in empirical food webs (i.e. trophic food chain, omnivory, exploitative and apparent composition)37,57 to assess a species motif census, although the general model would be applicable to any \(p\)-species motifs.
Network-scale cumulative effects scores (\({C}_{N}\)) were predicted in each grid cell as follows:
$${C}_{{N}_{x}}=\sum _{i}\frac{1}{\left|{M}_{i}\right|}\sum _{{m}_{i,x}\in {M}_{i}}\sum _{j}{D}_{j}\text{*}\overline{{\mu }_{j}}\text{*}{T}_{{i}_{{m}_{i,x}}}$$
where \(i\) is the focal species, \({M}_{i}\) is the motif census of species \(i\), \({m}_{i,x}\) are the 3-species motifs of interest forming species \(i\)’s motifs census, and \({D}_{j}\) is the log-transformed and scaled intensity of stressor \(j\).
\(\overline{{\mu }_{j}}\) corresponds to the joint sensitivity to stressor \(j\) of the species involved in motif \({m}_{i,x}\). Here, we explicitly consider that a species’ response to stressors depends on its own response as well as the response of species it interacts with. The joint sensitivity is measured as:
$$\overline{{\mu }_{j}}={w}_{1}{\mu }_{i,j}+{w}_{2}\sum _{k}^{2}{\mu }_{k,j}$$
\({\mu }_{i,j}\) and \({\mu }_{k,j}\) are the sensitivities to stressor \(j\) of focal species \(i\) and of the two species interacting with focal species \(i\) in motif \({m}_{i,x}\), respectively. \({w}_{1}\) and \({w}_{2}\) are weighting factors that give a relative importance to direct – i.e. effects to species \(i\) – and indirect – i.e. effects propagating through species \(k\) to species \(i\) – effects in the assessment. \({w}_{1}+2{w}_{2}=1\) to directly relate the weighting to a percent contribution to direct and indirect effects. For this assessment, we used \(w1=0.5\) and \(w2=0.25\) to give equal weight to direct and indirect effects.
\({T}_{{i}_{{m}_{i,x}}}\) is the trophic sensitivity of species \(i\) in motif \({m}_{i,x}\). It captures a species sensitivity to trophically-mediated effects, which depends on the structure of the community, the trophic position of focal species \(i\) and the specific entry points of stressors in the system14.
Stressors
We used the spatial distribution and intensity of 18 stressors available through an open-knowledge platform called eDrivers32. Stressors are divided into 4 groups: land-based (i.e. inorganic pollution, organic pollution, nutrient input, coastal development, and direct human impact), climate (i.e. positive and negative bottom-water and surface-water temperature anomalies, ocean acidification, and hypoxia), fisheries (i.e. demersal destructive, demersal non-destructive high-bycatch, demersal non-destructive low-bycatch, pelagic high-bycatch, and pelagic low-bycatch) and marine traffic (i.e. shipping and marine pollution; Extended Data Table 1). Methods to characterize each stressor are described in 32. Stressors with non-normal frequency distributions were log-transformed to avoid underestimating intermediate stressor intensity values26. All stressors were scaled between 0 and 1 to obtain relative intensities and allow comparisons between stressors. For each stressor, the 99th quantile of intensity distribution was used as the upper bound for scaling to control for extreme values that may or may not be real observations.
Species distribution
Biotic data
We used data from 4 monitoring programs conducted by Fisheries and Oceans Canada (DFO) to obtain a list of taxa with observed occurrences in the St. Lawrence marine ecosystem (Extended Data Table 2). We included a list of 30 known whale and seal species in the St. Lawrence and used distribution ranges available from the IUCN Species Red List of Threatened Species for 24 of the 30 marine mammal species (Extended Data Table 2)58. We curated the list of taxa used for the analyses by grouping and removing taxa based on expert knowledge and bibliographic research. For example, species of the same genus and hard to distinguish were grouped; species that were identified as probable misidentifications or outliers were removed from the data. The curation process yielded 424,953 taxa occurrences and 434,851 taxa absences for 391 taxa between 2010 and 2015 (Extended Data Table 2, Supplementary Table 1). The curation process was documented and is available on GitHub (https://github.com/eBiotic/Biotic). All species scientific names were resolved using the taxize R package59,60.
Abiotic data
We used environmental data characterizing the bottom-water and surface-water salinity, temperature, oxygen, primary productivity, pH (surface) and aragonite (bottom) conditions in the St. Lawrence marine ecosystem. We also considered latitude, longitude and depth, for 13 environmental descriptors. The data was accessed through various regional61–65 and global66 environmental monitoring programs and public repositories (Extended Data Table 3).
Spatial distribution
We extrapolated and mapped the distribution of taxa in the St. Lawrence marine ecosystem using the Random Forest ensemble learner67. We used the default parameters proposed by the randomForest R package to classify the presence or absence of taxa: 500 trees and the number of variables in the random subset at each node set to the square root of the number of variables68. We only considered taxa with at least 50 observations, yielding a total of 169 taxa (Supplementary Table 1). Each taxon was modelled using all 13 environmental descriptors (Extended Data Table 3). We generated pseudo-absences for taxa without absences in the dataset (\(n=5\)) by randomly sampling the study area at least 5 km away from observed points; for these taxa, we generated the same number of pseudo-absences as observed occurrences69. We measured the performance of the models for each taxon using the sensitivity, specificity, accuracy and True Skilled Statistics (TSS; Supplementary Table 1)70. For each taxon, we predicted spatial distribution within the same 1 \(k{m}^{2}\) resolution grid used for stressors. Individual taxa distributions were then smoothed using bisquare kernel smoothing71 with a 5 km radius to avoid potentially granular distributions that would affect estimations of species co-occurrence. We used distribution ranges available from the IUCN Species Red List of Threatened Species for 24 species of marine mammals58. Our dataset includes distribution maps for 193 taxa in the St. Lawrence marine ecosystem. We assumed that phytoplankton and zooplankton species were present throughout the St. Lawrence since these taxa are missing from our dataset and are required to consider trophic dynamics properly.
Species-specific sensitivity
Traits data
We documented the body composition, the maximal body size, the type of marine environment in which species are found, the feeding mode, the mobility and the phylum (Supplementary Table 2) of all 391 taxa available in the biotic dataset. We extracted traits data from the World Register of Marine Species (WoRMS)72, FishBase73, SeaLifeBase74, the Encyclopedia of Life75 and the Global Biotic Interaction (GloBI) database76,77. We used the taxize59,60, worrms78 and rfishbase79 R packages to extract traits data. Any taxon for which traits were unavailable programmatically were searched manually on the WoRMS and Encyclopedia of Life web portals. We also documented whether a species was targeted by fisheries or caught as bycatch by local fisheries using data from DFO’s Fisheries Logbook Program80.
Species-specific sensitivity
We evaluated the species-specific sensitivity of all 391 taxa to each stressor using a trait-matching approach. For each stressor, we identified traits that were known or suspected to influence a species sensitivity to the effects of the stressor (Supplementary Table 3). For example, the feeding strategy of an organism affects its sensitivity to nutrient and metal loading81, whereas its body composition affects its sensitivity to ocean acidification48. Traits were categorized to reflect their relative contribution to the sensitivity of a taxa to the effects of a stressor. For example, suspension feeders are generally more affected by nutrients and metals than deposit feeders81, whereas calcifying organisms are more vulnerable to the effects of ocean acidification than non-calcifying organisms48. Traits were categorized by giving a weight between 0 and 1 that reflects their relative contribution to the sensitivity of a taxa to a stressor: a weight of 0 represents a trait rendering taxa insensitive to the effects of a stressor, whereas a weight of 1 represents a trait associated with the highest relative sensitivity of a taxon to the effects of the stressor. The maximal sensitivity weight was retained if a taxon had multiple traits in a single category (e.g. crawler and swimmer). This sensitivity assessment was informed by expert knowledge and bibliographic research. Trait-matching rules and relative sensitivity weights for each stressor are available in Extended Data Table 5. The relative sensitivity of each taxon to a stressor was then evaluated as the product of the relative sensitivity weight of all traits associated with taxa sensitivity to the effects of the stressor. For example, the relative sensitivity to ocean acidification was evaluated using environment, mobility, body composition and phylum traits (Extended Data Table 5). This process yielded a relative sensitivity assessment ranging between 0 and 1 for each taxa.
Metaweb
We predicted the metaweb of the St. Lawrence marine ecosystem, i.e. the network of biotic interactions, using a recommender approach30. Here, we provide a brief overview of the approach, but refer to 30 for more details. The approach consists of a series of logical steps that predict a candidate resource list for each taxon based on empirical data available and the similarity among consumers and resources. It uses the K-nearest neighbour algorithm (KNN)82 to predict pairwise interactions given taxonomic and dietary similarity between consumers and resources and is informed by a catalogue of empirically known biotic interactions worldwide30. The interactions catalogue was built using food web data12,83,84, predator-prey interactions85 and pairwise interactions from the GloBI database76,77. We limited the compendium to taxa found in marine and coastal ecosystems. Taxa similarity was evaluated from taxonomic classification and sets of consumers or resources using the Tanimoto similarity measure. A weight of 0.5 was given to taxonomy and consumers or resources to consider them simultaneously86. The taxonomy of all taxa considered was accessed and validated from WoRMS72 using the taxize package59,60. We included the main phytoplankton and zooplankton taxa found in the St. Lawrence marine ecosystem to predict the metaweb87–89; we then grouped predictions under phytoplankton or zooplankton. This yielded a total of 393 taxa (\(S\)), considering all 391 taxa identified through the biotic data and the addition of phytoplankton and zooplankton. We predicted a metaweb structured by 4880 links (\(L\)), a link density (\({L}_{moy}=L/S\)) of 12.42 and a connectance (\(C=L/{S}^{2}\)) of 0.03, which is within range of most reported food webs90.
Trophic sensitivity
We provide a brief overview of the approach used by 14 to evaluate a species’ sensitivity to multiple stressors given its trophic position, i.e its trophic sensitivity. The effects of multiple stressors on the dynamics of the most empirically abundant 3-species motifs – i.e. tri-trophic food chain, omnivory, exploitative competition and apparent competition – were simulated using Lotka-Volterra models91. The dynamics of a single species regulated by density-dependent growth was also simulated; this control was included to consider disconnected species in the metaweb, which may arise due to insufficient data or because a species consumes detritus, bacteria or particulate organic matter. It was also used for the species-scale cumulative effects assessment so that disconnected species had the same cumulative effects results in the network-scale and species-scale assessments.
Negative effects of stressors were simulated by modifying combinations of equilibria equation parameters of population resource growth, mortality, attack and conversion rates (i.e. up to 9 parameters and 511 distinct pathways of effect)14. Modifications to parameters simulate effects of stressors on ecological processes; these represent the pathways through which stressors directly and indirectly affect ecological communities. The set of all ecological processes affected by stressors across species combine to collectively affect a community and form a pathway of effect. For each 3-species motif, all possible pathways of effects were simulated; this resulted in 127 unique pathways of effect for tri-trophic food chain, exploitative competition and apparent competition motifs through 7 parameters, 511 pathways of effect for the omnivory motif through 9 parameters, and 1 pathway of effect for disconnected species through 1 parameter. For each pathway of effect, a species’s trophic sensitivity was defined as the difference in its equilibrium abundance before and after the permanent appearance of stressors in the system; this represents the net effect of stressors on species and integrates all direct and indirect effects propagating to a focal species14,92.
Here, due to the challenge of empirically attributing effects to specific ecological processes, we simplified pathways of effects and broadly considered effects to a species density rather than effects to specific ecological processes as in 14. We used trophic sensitivities across possible pathways of effect simulated in 14 as heuristics to assess a species’ trophic sensitivity to the effects of stressors given its position in 3-species motifs. We used the absolute values of simulated trophic sensitivities and considered that any effect to a species’ population dynamics, whether negative or positive, can propagate and disturb the dynamics of an ecological community. We then simplified pathways of effects as a function of their effects to the density of each species. For example, a pathway of effect targeting the mortality of a consumer was considered to affect the density of that consumer, whereas a pathway of effect targeting the attack rate of a consumer was considered to affect the density of both the consumer and the resource. We averaged trophic sensitivities to pathways of effects according to their contribution to the effects on species density, resulting in 8 possible pathways of effect for all 3-species motifs and 2 pathways of effect for the single species, each with a pathway of effect where no effects are observed. Finally, we centered trophic sensitivities so that disconnected species had a value of 1, resulting in a trophic sensitivities ranging from 0 to 6.17, with an average of 1.67 \(\pm\) 1.57 and a median of 1.31.
Assessment and spatial data representation
The open-source software R 4.2.3 was used for all analyses93 and the package rcea94 was used to perform the cumulative effects assessment. See Extended Data Table 4 for a list of all R packages used. All datasets are presented at a 1 \(k{m}^{2}\) resolution even though some source data had coarser resolutions (Extended Data Table 1). We resampled and reprojected data when necessary using nearest neighbour estimates, which preserves the values of the source data. By doing so, we assume that the coarser data are evenly distributed across finer-scale cells with which they overlap. We used the NAD83 / Quebec Lambert projection (EPSG: 32198), which is well suited to represent and preserve surface area within our study system.
Methods references
57. Camacho, J., Stouffer, D. B. & Amaral, L. A. N. Quantitative analysis of the local structure of food webs. Journal of Theoretical Biology 246, 260–268 (2007).
58. IUCN. The IUCN Red List of Threatened Species. Version 2020-1. https://www.iucnredlist.org. Downloaded on 2020-07-01. (2020).
59. Chamberlain, S. A. & Szöcs, E. Taxize: Taxonomic search and retrieval in R. F1000Research 2, (2013).
60. Chamberlain, S. et al. Taxize: Taxonomic information from around the web. (2020).
61. Dutil, J.-D., Proulx, S., Chouinard, P.-M. & Borcard, D. A hierarchical classification of the seabed based on physiographic and oceanographic features in the St. Lawrence. Can. Tech. Rep. Fish. Aquat. Sci. 2916: Vii + 72 p. vii + 72 p (2011).
62. Dutil, J.-D. et al. Coastal and epipelagic habitats of the estuary and Gulf of St. Lawrence. Can. Tech. Rep. Fish. Aquat. Sci. 3009: Ix + 87 pp. ix + 87 pp (2012).
63. Galbraith, P. S. et al. Physical Oceanographic Conditions in the Gulf of St. Lawrence during 2017. DFO Can. Sci. Advis. Sec. Res. Doc. 2018/050. V + 79 p. v + 79 p. (2018).
64. Starr, M. & Chassé, J. Distribution of omega aragonite in the Estuary and Gulf of St. Lawrence in eastern Canada. Department of Fisheries and Oceans. (2019).
65. Blais, M. et al. Chemical and Biological Oceanographic Conditions in the Estuary and Gulf of St. Lawrence during 2017. DFO Can. Sci. Advis. Sec. Res. Doc. 2019/009. Iv + 56 p. iv + 56 p. (2019).
66. Assis, J. et al. Bio-ORACLE v2.0: Extending marine data layers for bioclimatic modelling. Global Ecology and Biogeography 27, 277–284 (2018).
67. Breiman, L. Random forests. Machine learning 45, 5–32 (2001).
68. Liaw, A. & Wiener, M. Classification and regression by randomForest. R news 2, 18–22 (2002).
69. Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: How, where and how many? Methods in Ecology and Evolution 3, 327–338 (2012).
70. Allouche, O., Tsoar, A. & Kadmon, R. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology 43, 1223–1232 (2006).
71. Dos Santos, A., Sémécurbe, F. & Pramil, J. Btb: Beyond the border - kernel density estimation for urban geography. (2022).
72. WoRMS Editorial Board. World Register of Marine Species. Available from http://www.marinespecies.org at VLIZ. Accessed 2018-09-19. https://doi.org/10.14284/170. (2017) doi:10.14284/170.
73. Froese, R. & Pauly, D. FishBase. World Wide Web electronic publication. Www.fishbase.org. Version (12/2019). (2019).
74. Palomares, M. & Pauly, D. SeaLifeBase. World Wide Web electronic publication. Www.sealifebase.org, version (12/2019). (2019).
75. Encyclopedia of Life. Encyclopedia of Life. Available from http://eol.org. Accessed 2020-06-01. (2020).
76. Poelen, J. H., Simons, J. D. & Mungall, C. J. Global biotic interactions: An open infrastructure to share and analyze species-interaction datasets. Ecological Informatics 24, 148–159 (2014).
77. Poelen, J., Gosnell, S., Slyusarev, S. & Waters, H. Rglobi: R interface to global biotic interactions. (2022).
78. Chamberlain, S. Worrms: World register of marine species (WoRMS) client. (2020).
79. Boettiger, C., Lang, D. T. & Wainwright, P. C. Rfishbase: Exploring, manipulating and visualizing FishBase data from R. Journal of Fish Biology 81, 2030–2039 (2012).
80. DFO. Zonal Interchange File Format (ZIFF) data. A compilation of landing data from logbook data between 2010 and 2015. Gestion des données, Institut Maurice Lamontagne, Department of Fisheries and Oceans (DFO) Mont-Joli, Canada. (2016).
81. Ellis, J. I. et al. Multiple stressor effects on marine infauna: Responses of estuarine taxa and functional traits to sedimentation, nutrient and metal loading. Scientific Reports 7, 12013 (2017).
82. Murphy, K. P. Machine learning : A probabilistic perspective. (MIT Press, 2012).
83. Brose, U. et al. Body sizes of consumers and their resources. Ecology 86, 2545–2545 (2005).
84. University of Canberra. Food Web Database - University of CANBERRA. (2016).
85. Barnes, C. et al. Predator and prey body sizes in marine food webs. Ecology 89, 881–881 (2008).
86. Desjardins-Proulx, P., Poisot, T. & Gravel, D. Ecological interactions and the Netflix problem. (2016).
87. Morissette, L. et al. Data gathering and input parameters to construct ecosystem models for the northern Gulf of St. Lawrence(mid-1980 s). Can. Tech. Rep. Fish. Aquat. Sci./Rapp. Tech. Can. Sci. Halieut. Aquat. 100 (2003).
88. Savenkoff, C. et al. Input data and parameter estimates for ecosystem models of the southern Gulf of St. Lawrence (mid-1980s and mid-1990s). 105 (2004).
89. Savenkoff, C. Input data and parameter estimates for ecosystem models of the lower St. Lawrence Estuary (2008). 150 (2012).
90. Dunne, J. A., Williams, R. J. & Martinez, N. D. Network structure and biodiversity loss in food webs: Robustness increases with connectance. Ecology letters 5, 558–567 (2002).
91. Gellner, G. & McCann, K. S. Consistent role of weak and strong interactions in high- and low-diversity trophic food webs. Nature Communications 7, 11180 (2016).
92. Abrams, P. A., Menge, B. A., Mittelbach, G. G., Spiller, D. A. & Yodzis, P. The Role of Indirect Effects in Food Webs. in Food Webs: Integration of Patterns & Dynamics (eds. Polis, G. A. & Winemiller, K. O.) 371–395 (Springer US, 1996). doi:10.1007/978-1-4615-7007-3_36.
93. Team, R. C. R: A Language and Environment for Statistical Computing. (2023).
94. Beauchesne, D. & Cazelles, K. Rcea: An R package to perform spatially-explicit cumulative effects assessments. (2023).
95. Beauchesne, D. eDrivers: Data on drivers of environmental change in the st. Lawrence system in eastern canada. (2020).
96. Ross, N. Fasterize: Fast polygon to raster conversion. (2022).
97. Bache, S. M. & Wickham, H. Magrittr: A forward-pipe operator for R. (2022).
98. Cazelles, K. & Beauchesne, D. Motifcensus: Count 3-nodes motifs for unipartite networks. (2023).
99. Hijmans, R. J. Raster: Geographic data analysis and modeling. (2023).
100. Xie, Y., Allaire, J. J. & Grolemund, G. R markdown: The definitive guide. (Chapman and Hall/CRC, 2018).
101. Xie, Y., Dervieux, C. & Riederer, E. R markdown cookbook. (Chapman and Hall/CRC, 2020).
102. Allaire, J. et al. Rmarkdown: Dynamic documents for R. (2023).
103. Pebesma, E. Simple features for r: Standardized support for spatial vector data. The R Journal 10, 439–446 (2018).
104. Pebesma, E. & Bivand, R. Spatial data science: With applications in R. 352 (Chapman and Hall/CRC, 2023).
105. Bivand, R. S., Pebesma, E. & Gomez-Rubio, V. Applied spatial data analysis with R, Second edition. (Springer, NY, 2013).
106. Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. R News 5, 9–13 (2005).
107. Pebesma, E. Stars: Spatiotemporal arrays, raster and vector data cubes. (2022).
108. Wickham, H. et al. Welcome to the tidyverse. Journal of Open Source Software 4, 1686 (2019).