# Integrating food webs in species distribution models can improve ecological niche estimation and predictions

Giovanni Poggiato, Jérémy Andréoletti, Laura J Pollock, Wilfried Thuiller

<figure><img src="image_1.png" /></figure>

## To cite this version:

Giovanni Poggiato, Jérémy Andréoletti, Laura J Pollock, Wilfried Thuiller. Integrating food webs in species distribution models can improve ecological niche estimation and predictions. Ecography, 2024, 10.1111/ecog.07546. hal-04985580

HAL Id: hal-04985580

https://hal.science/hal-04985580v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ECOGRAPHY

Research article

# Integrating food webs in species distribution models can improve ecological niche estimation and predictions

Giovanni Poggiato<fnref n="1" />,<fnref n="2" />, Jérémy Andréoletti<fnref n="1" />, Laura J. Pollock<fnref n="3" /> and Wilfried Thuiller<fnref n="1" />

**Correspondence:** Wilfried Thuiller (wilfried.thuiller@univ-grenoble-alpes.fr)

Ecography

2025: e07546

doi: 10.1111/ecog.07546

Subject Editor:
Catherine H. Graham
Editor-in-Chief: Miguel Araújo
Accepted 18 November 2024

<figure><img src="image_4.png" /></figure>

Biotic interactions play a fundamental role in shaping multitrophic species communities, yet incorporating these interactions into species distribution models (SDMs) remains challenging. With the growing availability of species interaction networks, it is now feasible to integrate these interactions into SDMs for more comprehensive predictions. Here, we propose a novel framework that combines trophic interaction networks with Bayesian structural equation models, enabling each species to be modeled based on its interactions with predators or prey alongside environmental factors. This framework addresses issues of multicollinearity and error propagation, making it possible to predict species distributions in unobserved locations or under future environmental conditions, even when prey or predator distributions are unknown. We tested and validated our framework on realistic simulated communities spanning different theoretical models and ecological setups. scenarios. Our approach significantly improved the estimation of both potential and realized niches compared to single SDMs, with mean performance gains of 8% and 6%, respectively. These improvements were especially notable for species strongly regulated by biotic factors, thereby enhancing model predictive accuracy. Our framework supports integration with various SDM extensions, such as occupancy and integrated models, offering flexibility and adaptability for future developments. While not a universal solution that consistently outperforms single SDMs, our approach provides a valuable new tool for modeling multitrophic community distributions when biotic interactions are known or assumed.

**Keywords:** Biodiversity modeling, biotic interactions, metanetworks, species distribution models, structural equation models

Introduction

Multitrophic interactions are generally known to shape species distributions, as illustrated by classic examples such as Gause's *Paramecium* experiments (Gause 1934) and the hare-lynx dynamics (MacLulich 1936). Yet, an important area of research is to

© 2025 The Author(s). Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_5.png" /></figure>

<footnote n="1">¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Grenoble, France</footnote>

<footnote n="2">²Univ. Grenoble Alpes, CNRS, INRIA, LJK, F-Grenoble, France</footnote>

<footnote n="3">³Department of Biology, McGill University, Montréal, Canada</footnote>
understand how strong and spatially variable are the effects
of multitrophic interactions on species (Thuiller et al. 2015),
and whether or not, there is a need for biodiversity models to
incorporate these effects for effective conservation planning
(Rosa et al. 2017, Pollock et al. 2020, Kass et al. 2024).

Species distribution models (SDMs, Guisan and Thuiller
2005) are widely used in biodiversity modeling, but they typ-
ically overlook biotic interactions. Over the past two decades,
considerable effort has focused on integrating biotic inter-
actions into models to improve predictions of species distri-
butions and biodiversity across spatial and temporal scales
(Guisan and Rahbek 2011, Warton et al. 2015a).

Emerging approaches, such as joint species distribution
models (JSDMs, Pollock et al. 2014, Clark et al. 2014,
Warton et al. 2015b) and Bayesian networks (Larsen et al.
2012, Ramazi et al. 2021), have been proposed to infer spe-
cies interactions from co-occurrence data. Yet, While JSDMs
offer statistical benefits for community-level data, they do
not explicitly infer interactions and often reflect correlations
arising from various unobserved factors (Zurell et al. 2018,
Poggiato et al. 2021). Bayesian networks mitigate some of
these issues but still rely on indirect inferences of species asso-
ciations from co-occurrence data, which has been criticized
for its limitations (Blanchet et al. 2020).

A more direct approach is to incorporate known or hypoth-
esized biotic interactions into SDMs. When documented,
these interactions can be added as supplementary covariates,
such as the presence of prey or competitors (Araújo and Luoto
2007, Thuiller et al. 2018). Studies have shown that incor-
porating biotic interactions can enhance model performance
and reveal how such interactions constrain species distribu-
tions, potentially distinguishing realized niches from funda-
mental niches (Gotelli et al. 2010, Boulangeat et al. 2012).
Historically, interactions were often assumed in a predefined
direction, which made this approach particularly useful for
modeling trophic interactions and, in some cases, competi-
tive exclusion (Leathwick and Austin 2001). However, this
method has rarely been extended to complex multitrophic
networks due to data limitations (i.e. the ‘Eltonian shortfall’,
Hortal et al. 2015) and technical complexities in predicting
interdependent species in multitrophic networks (Wisz et al.
2013). A preliminary attempt to account for known multi-
trophic interactions was made by Staniczenko et al. (2017,
2018), who used a Bayesian network to adjust SDM pre-
dictions based on predefined interaction networks. While
promising, this method required a fully known probability
table for interactions, making it challenging to apply to large
species networks. Although the authors proposed an optimi-
zation method to estimate the probability table, the compu-
tational demands limit its broader application. Additionally,
this approach views biotic interactions as a hierarchical filter
applied after environmental influences, rather than modeling
these factors jointly.

Models accommodating biotic interactions (Staniczenko et al. 2017) gain relevance with advances in interaction measurement (e.g. camera traps, gut content), open databases (e.g. GLOBI), and tools to query the literature

(Le Guillarme and Thuiller 2022, Grenié et al. 2023). These
growing datasets have accelerated the availability of large spe-
cies interaction networks, named metawebs (Dunne 2006).
These metawebs extend the concept of regional species pools
by incorporating interspecies interactions across trophic
levels (Albouy et al. 2019, Maiorano et al. 2020, Calderón-
Sanou et al. 2021, Potapov 2022). Gaps in knowledge can now
be addressed with trait-matching models (Pichler et al. 2020,
Caron et al. 2022) or phylogenetic analyses (Strydom et al.
2021, 2022). As we close the Eltonian shortfall, we can more
accurately use trophic interactions to develop ecologically
realistic predictive models (Windsor et al. 2022).

In this study, we introduce a flexible framework that
directly integrates metaweb information into SDMs using
structural equation modeling. Our framework models each
species based on its prey (for bottom-up control) or predators
(for top-down control) along with environmental variables,
allowing the prediction of community composition even at
sites where prey distributions are unknown. Implemented
in a Bayesian framework, our approach addresses multicol-
linearity and propagates uncertainty, enhancing predictions
for species strongly influenced by biotic factors. Additionally,
this framework can be used to test of ecological hypotheses
about biotic interactions’ roles in shaping species distribu-
tions and can identify key species exerting strong biotic
control.

We begin by describing our framework, discussing its
connections to existing models and outlining potential
extensions. We then evaluate its performance on simulated
ecological communities under a range of scenarios to identify
when it outperforms single SDMs. To promote accessibility
and encourage future developments, we have implemented
our framework in the R package ‘webSDM’, available on
CRAN (https://cran.r-project.org/web/packages/webSDM/
index.html).

Model framework

We present a method for modeling species distributions by
integrating both environmental variables and known spe-
cies trophic interactions. This approach first fits independent
models for each species (Fig. 1a–b), then combines them to
sequentially predict distributions for entire species pools,
including at unobserved sites where information on prey dis-
tributions may be missing (Fig. 1c).

To ensure that the models avoid circular dependen-
cies (which could lead to simultaneity bias; Pearl 2009),
the metaweb must be structured as a directed acyclic graph
(DAG). This requires selecting a directional control (e.g. bot-
tom-up or top-down), in which species are modeled either
by their prey distributions (bottom-up control) or by their
predators (top-down control). In this framework, we assume
1) a fully known, loop-free metaweb, 2) a stable metaweb
across space and time, and 3) a single directional control
specified in advance (though real ecosystems may exhibit
both types of control).
<figure><img src="image_6.png" /><figcaption>Figure 1. Description of the proposed model. The input data are the environmental variables, species distribution data, and a directed acyclic graph representing the species trophic interactions (a). Assuming six hypothetical species, each species is modeled as a function of the environmental covariates and its prey (or predators). Since the likelihood factorizes, each species can be modelled independently (b). Predictions under new environmental conditions (in space or time) where prey distribution are unknown (c). To predict with our trophic SDM, we need to guarantee that, for each predator, all its prey species have already been predicted. To do so, we order species according to the topological order σ derived from the metaweb and then perform predictions following this order. Therefore, predictions of species Y₄, Y₅ and Y₆ are obtained as a function of the environmental conditions x_new. Then, predictions for Y₂ are computed as a function of x_new and the predictions of Y₄ and so on for species Y₁ and Y₅.</figcaption></figure>

## Model fitting

Given a metaweb structured as a directed acyclic graph G and an occurrence matrix Y (where each row i represents a sampling unit and each column j a species), we model the distribution of each species as a function of environmental covariates and its prey. The environmental covariate matrix X has columns indexed by k = 1, ..., p. For each species j, we use a generalized linear model (GLM), defined as follows:

$$ 
\begin{gathered}
y_{ij} \sim f(\mu_{ij}) \\
g(\mu_{ij}) = \beta_{j0} + \sum_{k=1}^{p} \beta_{jk} x_{ik} + \sum_{j':(j',j) \in E} \alpha_{jj'} y_{ij'} \quad (1)
\end{gathered}
$$

where f() is the probability distribution of the observed species (e.g. binomial for presence-absence or Poisson for counts) with parameter μᵢⱼ (we omitted the eventual dispersion parameter in the equation for simplicity) and g() is the corresponding link function. βⱼₖ is the effect of environmental covariate k on species j (with βⱼ₀ as the intercept of species j), and αⱼⱼ represents the influence of prey j' on predator j, and E is the set of links of the metaweb G. This structure ensures conditional independence within the model, where each species depends on its immediate prey rather than on prey further down the trophic chain (Pearl's 'causal Markov condition'; Pearl 2009). Consequently, the joint likelihood factorizes:

$$ p(\mathbf{y}|\mathbf{x}) = \prod_{j=1}^{S} p(\mathbf{y}_j | \mathbf{X}, \text{prey}(\mathbf{y}_j)) \quad (2) $$

This is not a model assumption, but a mathematical property of the model arising from Eq. 1. It implies that we can estimate the whole model as a set of independent GLMs, one for each species (Grace et al. 2012). Interestingly, modeling species as a function of their prey (or their predators) and the environmental covariates is equivalent to a structural equation model (Shipley 2000, Grace 2006) whose structure is given by G, plus a dependence of all species to the environment (cf. Box 1 for a description of SDM, JSDM, Bayesian networks, and our trophic model through the lens of SEMs). The estimation of the parameters of a SEM (whose structure is a directed acyclic graph) as a set of independent GLMs is a common practice and takes the name of 'local estimation' or piecewise SEM (Shipley 2000, Lefcheck 2016). Therefore, model fitting simply requires fitting S local models, here GLMs, which can be computationally fast and can even be parallelized. Thus, computational requirements grow linearly with the number of species (without parallelization), while the effect of the number of sites depends on the chosen technique to fit GLMs.

If the model described in Eq. 1 can be fitted separately for each species, the ensemble of species distribution models has then to be combined to generate predictions on the whole species pool at unobserved sites, given the metaweb.

## Model predictions

Predicting across species pools using SEMs can be complex, as predictions for higher trophic levels require prior knowledge of prey distributions. We address this using the metaweb's topological ordering, a sequence that ensures all prey are predicted before their respective predators. By predicting in this
**Box 1. Structural dependence in species distribution modeling**

Like any multivariate model, all kinds of species distribution models introduce probabilistic dependencies between vari-
ables (i.e. both species and environmental covariates). These dependencies can be expressed using the graphical representa-
tion of structural equation models (SEM), i.e. the ‘path diagram’, to understand the differences between models and, in
particular, the relationships they introduce between species. In SEM, each variable is represented as a box, a direct arrow
from one box to another indicates a direct relationship between the two variables (e.g. an arrow from A to B means that
A affects B), and an undirected arrow between two variables indicates a correlation. From such a graphical representation,
we can depict the directed and indirect relationships between variables, as well as the marginal and conditional dependen-
cies that these relationships imply (cf. Koller and Friedman 2009 for a thorough introduction to probabilistic graphical
models). Here, we express the dependencies introduced by the single SDM, JSDM, Bayesian networks, and trophic SDM
when modeling a hypothetical set of six interacting species (Y₁ to Y₆) and a set of environmental covariates X (Box Fig. 1).

SDM

Single SDMs model each species independently as a function of the environment (Guisan and Thuiller 2005). This is therefore equivalent to an SEM diagram with arrows pointing from environment X to all species (Box Fig. 1a). As a consequence, species are marginally dependent, as they are correlated through their response to environment (e.g. species with similar niches will be predicted to co-occur), but they are conditionally independent given the environment. Since the effect of interacting species is not controlled for, the species-environment relationships capture both the environmental and biotic effect, and, as such, SDMs only infer species’ realized niches (i.e. biotic interactions are implicitly taken into account, Araujo and Guisan 2006).

JSDM

JSDMs model each species as a function of the environment assuming residuals are correlated across species (Pollock et al.
2014, Ovaskainen et al. 2017). This corresponds to an SEM diagram where the environment X points to all species, and

<figure><img src="image_7.png" /></figure>

Box Figure 1. Structural dependence of SDM (a), JSDM (b), and trophic SDM (c). Black boxes refer to the environment covariates (X), while black circles represent species (Y). Gray circles represent residuals (ε). Bayesian network correspond to the same diagram of trophic SDM (c), except that for Bayesian networks the graph is inferred from the data with obvious uncertainties.
species are linked with each other, through their residuals, with undirected arrows representing their correlations (Box Fig. 1b). These correlations imply a symmetrical, undirected, relationship between species residuals, so that species are conditionally dependent given the environment. Moreover, since conditional independencies are expressed in a latent layer, every pair of species is conditionally dependent given the rest of the network (i.e. the inferred graph at the species level is fully connected). However, due to the residual nature of these correlations, they have little effect on the estimates of the species-environment relationships (Chib 1998, Poggiato et al. 2021) and do not modify species marginal predictions (Poggiato et al. 2021).

## Bayesian network models

Bayesian network models (Larsen et al. 2012, Ramazi et al. 2021) infer a directed acyclic graph (DAG) from species distributions and environmental covariates. Once this network of co-occurrences is inferred, species are modeled as a function of their parents (similarly to the trophic SDM), depending on the method used (Montesinos-Navarro et al. 2018, Ramazi et al. 2021). The SEM diagram of a Bayesian network is then simply the inferred DAG. The advantage of Bayesian networks over JSDM is that they model the species together with the environment, so the two effects can be properly separated. However, the inferred DAG does not necessarily correspond to the true interaction networks, which leads to a conceptual difference between the Bayesian network and trophic SDM. Moreover, to our knowledge, Bayesian networks have never been used to generate species predictions at unobserved sites where prey distribution is unknown.

## Trophic SDM

Our approach proposes to inject the knowledge of the metaweb by modeling species as a function of their prey (or predators) and of the environment. This introduces a direct arrow from prey to predators on top of the arrows from the environment to species (Box Fig. 1c). So, the effect of prey is controlled for and the distribution of prey directly determines the prediction of the predators (vice versa for a top-down control). Each species is therefore conditionally independent, given its prey and the environment, to the preys of its preys. However, for any connected metaweb, species are marginally correlated given the environment only, so that, for example, two predators feeding on the same prey are marginally correlated due to the indirect effect of sharing the same prey (i.e. they tend to co-occur because they both feed on the same prey).

order, we create a 'cascade' from basal species up to top predators, capturing the interdependencies in predictions (Fig. 1c). This approach can be reversed for top-down scenarios.

Implemented in a Bayesian framework, our model samples from the posterior predictive distribution of each species, propagating uncertainty across predictions and generating credible intervals that widen progressively up the trophic chain (this allows to correctly estimate species' predictive posterior distribution, Supporting information). While computational demands grow with species count in prediction steps, the precise scaling depends on the number of posterior samples and metaweb complexity. In the further section, we describe alternative pathways to reduce computational times.

As a side note, the Bayesian (versus frequentist) implementation of our GLMs should not be confused with the afore-mentioned Bayesian networks, which are probabilistic graphical models that describe conditional dependencies among probabilistic variables.

## Dimension reduction: sparse modeling and composite variables

Incorporating species as predictors inherently introduces more predictors into the models, potentially leading to multicollinearity issues, particularly in large networks and for generalist species. While multicollinearity may not directly compromise model predictive performance, it can distort coefficient estimates, affecting our understanding of

prey-predator effects (Dormann et al. 2013, Tredennick et al. 2021).

Including non-influential predictors may also lead to overfitting, a risk in complex ecosystems where biotic interactions may not significantly impact species distributions at broad spatial scales (Pottier et al. 2013, Thuiller et al. 2015).

Both multicollinearity and overfitting can be mitigated by constraining coefficients of unimportant predictors to zero, reducing model complexity – a technique known as sparse modeling or regularization (O'Hara and Sillanpää 2009, Hastie et al. 2015). We implemented this sparse approach using the horseshoe prior (Carvalho et al. 2010), a method that heavily shrinks coefficients of less influential predictors while preserving stronger signals.

This approach significantly reduces model complexity and should guarantee that in cases where species interactions leave no signals in the data, our model should be equivalent to a single SDM. In extreme multicollinearity cases, the sparse model may randomly select between collinear variables, so we recommend comparing predictions from single and trophic SDMs to ensure interpretability.

For some generalist predators, there might simply be an intractable number of prey. In these cases, assuming that every prey has a differential effect on predator distribution is not only a statistical problem but might seem ecologically unjustified. Instead, we could expect the richness or diversity of prey, or whether at least a prey is available, to be important. Hence, we implemented the use of composite variables
that summarize the information of a large number of vari-
ables from the graph into a few summary variables (Henseler
2021). Implemented examples include using prey richness,
a binary variable representing prey availability, or grouping
prey into trophic clusters based on functional similarities
(Gauzens et al. 2015, O'Connor et al. 2020). Composite
variables streamline predictor sets, facilitating model inter-
pretation for generalists and top predators.

Implementation

Our framework is implemented in the R package ‘webSDM’
available on CRAN (https://cran.r-project.org/web/packages
/webSDM/index.html). A comprehensive illustration of
webSDM is provided in the Supporting information. The
independent GLMs are fitted using the R library ‘rstanarm’
(Goodrich et al. 2022), which enables fast sampling using the
STAN platform and offers diverse prior settings. Currently,
webSDM supports presence–absence data and includes
options for regularization, composite variables and various
model evaluation metrics. The ability of our algorithm to
recover true parameters under model-based sampling is vali-
dated in the Supporting information.

Validation against realistic simulated data

To validate our model, we tested it on realistic simulated
datasets that reflect ecological processes without excessive

stochasticity or sampling biases, following established prac-
tices (Zurell et al. 2010, Meynard et al. 2019). Using a
theoretical model, we simulated species distributions for
interacting species across varied environments, capturing
both bottom-up and top-down controls (Fig. 2a–b). We
then applied single-species distribution models (SDMs) and
trophic SDMs to the simulated data to compare their pre-
dictive performance for realized niches (observed species dis-
tributions) and fundamental niches (expected distributions
when biotic constraints are removed). By varying parameters
of the theoretical model, such as metaweb complexity and
species niche breadth, we explored how well the models per-
form under different ecological scenarios (Fig. 2c–d).

Simulation settings

Species communities were simulated using a generalized
Lotka–Volterra model to capture realistic ecological interac-
tions. We also run alternative simulations using variations of
this model and the Ricker model.

The Lotka–Volterra model describes species abundances
over time with ordinary differential equations, where prey
positively influence predators and predators negatively affect
prey (Lotka 1920, Volterra 1926). We included abiotic con-
trol on species growth rates along a single environmental gra-
dient. Following niche theory, we assumed that each species’
growth rate decreases in a Gaussian fashion as it deviates from
its optimal position along the gradient (Hirzel and Le Lay
2008). Species niche optima were randomly assigned, and

<figure><img src="image_8.png" /><figcaption>Figure 2. Application of the virtual ecologist approach to trophic SDMs (Zurell et al. 2010). The virtual ecological model (a) consists of defining, for a given number of species, their abiotic preferences as well as a trophic interaction network, and then making their abundances vary according to a selected theoretical model combining abiotic (blue) and biotic (red) controls. After a certain time (i.e. at the stationary state), species abundances are transformed in presence-absences (b) and used to fit the SDM and trophic SDM (c). Finally, the models are evaluated (d) by looking at the congruence of predicted realized and fundamental niches (dashed lines, and continuous lines represent the true niches).</figcaption></figure>
interactions were fixed along the gradient as per a directed
metaweb structure, with prey and predator impacts modeled
as asymmetric effects. To enable comparison between true
and predicted occurrence probabilities from single and tro-
phic SDM, we introduced a stochastic component on the
growth rate of the Lotka–Volterra model (detailed in the
Supporting information). For a given directed acyclic graph
and niche configuration, we then generated multiple com-
munities along the environmental gradient from the sto-
chastic generalized Lotka–Volterra model. Simulated species
abundances were then transformed to presence–absence (a
species was set as presence if its continuous abundance was
greater than zero) to obtain a species distribution dataset.

To test robustness across simulation parameters, we varied the species pool size (10, 20 or 50 species), the link probability in the metaweb (0.2, 0.3 or 0.5), and species' niche breadths (scaled by species pool size, set at 2, 6 or 10). Additional technical parameters, such as the number of points along the environmental gradient and the number of communities simulated for each environment, were also varied. Using Latin hypercube sampling, we selected 50 combinations of these simulation parameters. Then, for a given set of parameters, we sampled one interaction graph and species niches and used the above-described procedure to obtain one dataset of species distribution. This procedure was repeated 100 times, each of them with a randomly sampled metaweb and species niches, for each of the 50 parameter combinations. We thus obtained a total of 100 × 50 = 5000 species distribution datasets.

Statistical modeling and data analysis

The statistical analysis of the simulated data was carried out
independently for each of the 5000 species distribution data-
sets (Supporting information). Using the true interaction
network as the metaweb, we fitted our trophic SDM with a
second-degree polynomial for the environmental effect. We
also fitted generalized linear models (GLMs) for each species
as single SDMs, using the same polynomial term, allowing
direct comparison with the trophic model by excluding biotic
interaction terms. This ensured any performance differences
were due to the inclusion of trophic interactions and not
algorithmic variation.

We evaluated model performances through cross-valida-
tion, randomly dividing the 51 points along the environmen-
tal axis into five folds (i.e. training on four folds and predicting
on the remaining fold, repeated five times). Predictive perfor-
mance was assessed for both realized and fundamental niches.
We estimated the species' realized niche as its modelled prob-
ability of presence along the environmental gradient (sensu
Hutchinson 1957), which we compared to the observed
probability of presence. We defined the fundamental niche
as the environmental conditions where a species could poten-
tially be present if the constraints imposed by biotic inter-
actions were released (Soberón 2007). To release the biotic
constraints, we expressed the species' fundamental niche as
the probability of presence along the environmental gradient

when its prey is assumed to be present and its predators absent
everywhere. With our trophic model, we can predict species'
fundamental niches by setting the prey as present everywhere
across the environmental gradient (i.e. a conditional predic-
tion). Single SDMs cannot condition on the value of other
species, so we considered that the predictions of the potential
and realized niche coincide. We then compared the predicted
fundamental niches to the theoretical simulated ones, com-
puted as the probability of species having a positive growth
rate given the presence of prey, and the absence of predators
(Supporting information).

We evaluated how well the models predicted the real-
ized niche by comparing the predicted probability of pres-
ence to observed presence-absences using the area under
the ROC curve (AUC, Fawcett 2006). We also used the
Wasserstein distance with p=1 (to compare the distance
between the predicted and observed probability of pres-
ence), an approximation of the leave-one-out cross-valida-
tion likelihood (Vehtari et al. 2017), and calibration (i.e.
the number of times the 95% credible interval correctly cov-
ers the true value, Norberg et al. 2019). To measure model
performances in retrieving the fundamental niche, we only
used the Wasserstein distance and calibration (Supporting
information).

Finally, we checked the width of the 95% credible interval
of the predicted realized niche, to understand if our uncer-
tainty propagation technique could lead to overly large cred-
ible intervals, especially for top predator species.

The codes to fully reproduce the results are available in the
Supporting information.

Results

Overall, the trophic SDM effectively captured species' real-
ized niches, achieving an average AUC of 0.85. Across all
simulation conditions, the trophic model improved AUC by
an average of 6% compared to the SDM (one-sided paired
t-test, *p* < 10<sup>-16</sup> for non-basal species, Fig. 3a). Specifically,
the trophic model enhanced AUC for 66% of species, while
for others, predictions were only slightly worsened (with just
0.3% of species seeing an AUC decrease over 20%). The larg-
est improvements occurred when SDMs struggled to predict
realized niches accurately (Supporting information).

Regression analysis of AUC relative improvement as a
function of simulation parameters showed that the size of
the species pool and the normalized niche breadth (a stan-
dardized measure of species niche breadth) were the main
drivers of variation in model improvement (Supporting
information). Larger species pools generally boosted model
performance, while narrower niche breadths slightly reduced
it, suggesting the trophic model performed especially well
for species with specialized environmental requirements. To
illustrate these patterns, we examined two example species
across two simulation scenarios with differing niche breadths
but identical parameters otherwise (Fig. 4). In cases of nar-
row niche breadths, prey distributions were restricted to a
<figure><img src="image_9.png" /><figcaption>Figure 3. Model performances for the realized niche (a) and fundamental niche (b). Different colors correspond to different statistical models (SDM in red, trophic SDM in blue). In (a), AUC measures model performances in predicting the realized niche (higher AUC corresponds to improved model performances). In (b), the Wasserstein distance measures model performances in predicting the fundamental niche (lower Wasserstein distance corresponds to improved model performances).</figcaption></figure>

narrow environmental range, which produced a marked shift
in realized niche shape. Here, the SDM performed poorly,
while the trophic model provided a much more accurate pre-
diction (Fig. 3a, AUC: SDM=0.68, trophic SDM=0.86).
Conversely, when prey were distributed across a broader envi-
ronmental range, predators' distributions were less impacted
by biotic interactions, leading SDMs to closely match the
trophic model predictions (Fig. 3b, AUC: SDM=0.91, tro-
phic SDM=0.90).

The trophic model also achieved better predictions for spe-
cies' fundamental niches, reducing the Wasserstein distance to
the true fundamental niche by 8% on average (one-sided paired
t-test, p < 10<sup>-16</sup> for non-basal species, Fig. 3b, Supporting
information). Unlike realized niche prediction, this improve-
ment in fundamental niche prediction was not affected by any
simulation parameters. The most substantial gains were again
observed in cases where SDMs failed to predict the funda-
mental niche accurately (Supporting information).

These findings were consistent across other evaluation
metrics, including alternative theoretical models and param-
eter combinations. Detailed results for each evaluation metric
are available in the Supporting information.

As expected, the width of the credible intervals for the pre-
dicted probability of presence was greater for non-basal spe-
cies, increasing with prey count. However, this enlargement
stabilized around 0.25, regardless of trophic level and prey
count, indicating that credible interval widths remained prac-
tical across all network conditions (Supporting information).

Limitations and extension of the simulation
framework

Our validation tested a wide range of simulation parameters to establish the robustness of our method under idealized conditions. However, extending the simulations to capture real-world limitations would provide valuable insights. One natural extension would involve testing multiple environmental gradients instead of a single one to enhance ecological relevance. Additionally, exploring the effects of inaccuracies in the metaweb would quantify how error rates in network structure impact model performance.

Discussion

In this study, we introduced a flexible species distribution
modeling (SDM) framework that accounts for known trophic
interactions. Our framework models species distributions as
functions of both environmental variables and trophic rela-
tionships (e.g. prey or predator presence). It is designed to
handle multicollinearity and propagate error, enabling pre-
dictions for species at unobserved sites or under future con-
ditions where prey or predator distributions are unknown.
We rigorously validated this model with simulations that
included simultaneous top-down and bottom-up controls,
a realistic condition in many ecosystems. Despite these com-
plexities and the dynamic ecological properties simulated,
<figure><img src="image_10.png" /><figcaption>Figure 4. True and predicted realized niches of two species corresponding to simulations where species had low normalized niche breadth (a) and high normalized niche breadth (b). For species (a), trophic SDM improved predictions (AUC=0.86 for trophic SDM and 0.68 for SDM). For species (b) performances were almost equal (AUC=0.91 for SDM, 0.90 for trophic SDM).</figcaption></figure>

our model improved on single SDMs for predicting species'
realized niches, particularly when biotic interactions strongly
influenced these niches, and it also allowed for better estima-
tion of species' fundamental niches.

However, this framework is not universally superior to
SDMs, and it may not consistently outperform them. The
influence of biotic interactions on species distributions is
still debated, and potential inaccuracies in data sources (e.g.
metawebs, environmental covariates) could mean that includ-
ing prey as additional predictors may not always enhance
predictions and may sometimes degrade them. While regular-
ization helps manage these issues, making trophic SDM pre-
dictions at least equivalent to single SDMs, practical results
can vary. We therefore recommend using both SDMs and
trophic SDMs and selecting the better-performing model in
each case. Our sensitivity analysis highlighted conditions that
may enhance trophic SDM performance, but further testing
in real-world scenarios is necessary to refine these insights.

Model assumptions and comparison with (J)SDM

A primary assumption of our model is that the metaweb of
interactions remains static across environmental gradients,
both spatially and temporally. However, species often switch
prey when availability changes, as seen in cases of biological
invasions that involve new predators or prey. To address this
potential limitation, it is important to construct metawebs
that include all possible interactions, not just those observed
in specific contexts, as is common practice in published
metawebs. (Maiorano et al. 2020, Calderón-Sanou et al.
2021).

Modeling biotic interactions is inherently challenging
(Thuiller et al. 2013), and while causal relationships are

difficult to determine with static approaches (Blanchet et al.
2020), we have shown that explicitly incorporating known
interactions can improve predictive performance for realized
niches and even enhance our understanding of fundamental
niches. In the conditional prediction stage, the environment
ultimately determines species presence through its influence
on interacting species, meaning species respond indirectly to
the environment via the cumulative responses of their prey
(cf. the Supporting information for further discussion). In
contrast, highly flexible SDM algorithms (e.g. XGBoost,
RandomForest) with many environmental predictors may
implicitly capture interactions but risk overfitting, potentially
compromising predictions in unobserved areas. Our trophic
model, by comparison, allows for complex environmental
responses within a more parsimonious and ecologically inter-
pretable framework, especially when trained on observed data
(Supporting information).

Beyond a predictive model

Our framework is also a tool for investigating the role of
trophic interactions in species distributions and community
composition. It allows exploration of how trophic interac-
tion effects vary with spatial resolution (Thuiller et al. 2015),
trophic level, generalism (Fraser et al. 2021), functional traits
(Boet et al. 2020, van der Merwe et al. 2021), and environ-
mental conditions. Including interactions between abiotic
and biotic terms (Paquette and Hargreaves 2021) enables fur-
ther insights into environmental impacts on biotic relation-
ships. While we have focused on bottom–up interactions, a
complementary top–down analysis could clarify the direction
and strength of trophic control across systems. Additionally,
exploring various model specifications of biotic terms, such
as incorporating prey diversity, richness, or abundance, could reveal new facets of species interactions and control mechanisms.

This framework can also be valuable in studying network robustness under global change, predicting how species extinctions may cascade through the trophic chain. In this way, it provides a method to assess the potential impacts of climate and land use change on species distributions, community composition and trophic diversity.

## Model extensions

Our current framework is limited to presence-absence data, but species presence often depends not just on the presence of prey but also on prey abundance. Improvements in model accuracy may be possible by modeling abundance instead, for instance through a Poisson or negative binomial link function in generalized linear models (GLMs). This modification is straightforward and could yield more accurate results.

The trophic SDM is implemented in a Bayesian framework, which allows for error propagation and the integration of informative priors based on ecosystem knowledge (e.g. interaction network uncertainties). Alternatively, a frequentist approach could offer faster computations and compatibility with other statistical techniques. Local model estimation enables integration of diverse SDM extensions, such as presence-only data modeling (Renner et al. 2015), integration of multiple data sources (Isaac et al. 2020), or adjustments for imperfect species detection (MacKenzie et al. 2006). Though not yet implemented in our package, these extensions would be simple to incorporate given the local estimation approach. However, a drawback of local estimation is the inability to model species-environment coefficients hierarchically, which can be useful for modeling rare species (Pollock et al. 2014). A potential solution would be to adopt a global estimation approach, though this assumes all species respond to a uniform set of environmental variables – a questionable assumption for large ecological networks or diverse species with differing ecological traits (e.g. parasites, autotrophs, herbivores).

## Conclusions and perspectives

We believe our framework provides a promising approach to incorporating trophic networks into species distribution models. Currently, the model is limited to directed acyclic trophic networks, preventing it from handling competitive interactions, a key mechanism in trophic communities. However, there are ways to address this. For example, asymmetrical networks can be created by selecting a dominant interaction direction, allowing the framework to extend to non-trophic interactions such as competition or mutualism (e.g. plant-pollinator and plant-herbivore interactions; see Thuiller et al. 2018 for a single herbivore case). For cyclical or symmetrical networks, a promising approach involves

averaging predictions across the set of spanning trees in the network (Chaiken and Kleitman 1978), an ensemble technique already used in machine learning (Read et al. 2021). This method would enable integration of interaction uncertainties within the metaweb, though it requires significant development. We hope this modeling framework will stimulate further research and inspire advancements in the integration of ecological networks into species distribution models, pushing the boundaries of predictive ecology.

**Acknowledgements** – We thank Marc Ohlmann and Stéphane Robin for the insightful discussions around network topology and graphical model, and Julyan Arbel for his help with the estimation of the posterior predictive distribution.

**Funding** – This study has received funding from the ERA-Net BiodivERsA-Belmont Forum, with the national funder Agence Nationale de la Recherche (ANR-18-EBI4-0009), part of the 2018 Joint call BiodivERsA-Belmont Forum call (project ‘FutureWeb’). WT and GP also acknowledge support from the European Union’s Horizon Europe under grant agreement number 101060429 (project NaturaConnect), and the French Agence Nationale de la Recherche (ANR) through the Gambas (ANR-18-CE02-0025), EcoNet (ANR-18-CE02-0010), MIAI@Grenoble (ANR-19-P3IA-0003) projects.

## Author contributions

**Giovanni Poggiato:** Conceptualization (equal), Data curation (lead), Formal analysis (lead), Investigation (lead), Methodology (lead), Software (lead), Validation (lead), Visualization (lead), Writing - original draft (lead), Writing - review and editing (lead). **Jérémy Andréoletti:** Methodology (supporting). **Laura J. Pollock:** Conceptualization (supporting), Methodology (supporting), Supervision (supporting), Writing - original draft (supporting), Writing - review and editing (supporting). **Wilfried Thuiller:** Conceptualization (equal), Funding acquisition (lead), Investigation (equal), Methodology (equal), Project administration (lead), Resources (lead), Supervision (lead), Validation (equal), Writing - original draft (equal), Writing - review and editing (lead).

## Transparent peer review

The peer review history for this article is available at https://www.webofscience.com/api/gateway/wos/peer-review/ecog.07546.

## Data availability statement

Data are available from the Github Repository: https://github.com/giopogg/webSDM/tree/main/publication (Poggiato et al. 2024).

## Supporting information

The Supporting information associated with this article is available with the online version.
## References

Albouy, C., Archambault, P., Appeltans, W., Araújo, M. B., Beauchesne, D., Cazelles, K., Cirtwill, A. R., Fortin, M.-J., Galiana, N., Leroux, S. J., Pellissier, L., Poisot, T., Stouffer, D. B., Wood, S. A. and Gravel, D. 2019. The marine fish food web is globally connected. – Nat. Ecol. Evol. 3: 1153–1161.

Araújo, M. B. and Guisan, A. 2006. Five (or so) challenges for species distribution modelling. – J. Biogeogr. 33: 1677–1688.

Araújo, M. B. and Luoto, M. 2007. The importance of biotic interactions for modelling species distributions under climate change. – Global Ecol. Biogeogr. 16: 743–753.

Blanchet, F. G., Cazelles, K. and Gravel, D. 2020. Co-occurrence is not evidence of ecological interactions. – Ecol. Lett. 23: 1050–1063.

Boet, O., Arnan, X. and Retana, J. 2020. The role of environmental vs biotic filtering in the structure of European ant communities: a matter of trait type and spatial scale. – PLoS One 15: e0228625.

Boulangeat, I., Gravel, D. and Thuiller, W. 2012. Accounting for dispersal and biotic interactions to disentangle the drivers of species distributions and their abundances. – Ecol. Lett. 15: 584–593.

Calderón-Sanou, I., Münkemüller, T., Zinger, L., Schimann, H., Yoccoz, N. G., Gielly, L., Foulquier, A., Hedde, M., Ohlmann, M., Roy, M., Si-Moussi, S. and Thuiller, W. 2021. Cascading effects of moth outbreaks on subarctic soil food webs. – Sci. Rep. 11: 15054.

Caron, D., Maiorano, L., Thuiller, W. and Pollock, L. J. 2022. Addressing the Eltonian shortfall with trait-based interaction models. – Ecol. Lett. 25: 889–899.

Carvalho, C. M., Polson, N. G. and Scott, J. G. 2010. The horse-shoe estimator for sparse signals. – Biometrika 97: 465–480.

Chaiken, S. and Kleitman, D. J. 1978. Matrix tree theorems. – J. Comb. Theor. A 24: 377–381.

Chib, S. 1998. Analysis of multivariate probit models. – Biometrika 85: 347–361.

Clark, J. S., Gelfand, A. E., Woodall, C. W. and Zhu, K. 2014. More than the sum of the parts: forest climate response from joint species distribution models. – Ecol. Appl. 24: 990–999.

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D. and Lautenbach, S. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. – Ecography 36: 27–46.

Dunne, J. A. 2006. The network structure of food webs. – In: Pascual, M. and Dunne, J. A. (eds), Ecological networks: linking structure and dynamics. Oxford Univ. Press, pp. 27–86.

Fawcett, T. 2006. An introduction to ROC analysis. – Pattern Recognition. Lett. 27: 861–874.

Fraser, D., Soul, L. C., Tóth, A. B., Balk, M. A., Eronen, J. T., Pineda-Munoz, S., Shupinski, A. B., Villaseñor, A., Barr, W. A., Behrensmeyer, A. K., Du, A., Faith, J. T., Gotelli, N. J., Graves, G. R., Jukar, A. M., Looy, C. V., Miller, J. H., Potts, R. and Lyons, S. K. 2021. Investigating biotic interactions in deep time. – Trends Ecol. Evol. 36: 61–75.

Gause, G. F. 1934. Experimental analysis of Vito volterra's mathematical theory of the struggle for existence. – Science 79: 16–17.

Gauzens, B., Thébault, E., Lacroix, G. and Legendre, S. 2015. Trophic groups and modules: two levels of group detection in food webs. – J. R. Soc. Interface 12: 20141176.

Goodrich, B., Gabry, J., Ali, I. and Brilleman, S. 2022. rstanarm: Bayesian applied regression modeling via Stan. – http://mcstan.org/rstanarm/.

Gotelli, N. J., Graves, G. R. and Rahbek, C. 2010. Macroecological signals of species interactions in the Danish avifauna. – Proc. Natl Acad. Sci. UA 107: 5030–5035.

Grace, J. B. 2006. Structural equation modeling and natural systems. – Cambridge Univ. Press.

Grace, J. B., Schoolmaster, D. R., Jr, Guntenspergen, G. R., Little, A. M., Mitchell, B. R., Miller, K. M. and Schweiger, E. W. 2012. Guidelines for a graph-theoretic implementation of structural equation modeling. – Ecosphere 3: art73.

Grenié, M., Berti, E., Carvajal-Quintero, J., Dädlow, G. M., Sagouis, A. and Winter, M. 2023. Harmonizing taxon names in biodiversity data: a review of tools, databases and best practices. – Methods Ecol. Evol. 14: 12–25.

Guerra, C. A., Delgado-Baquerizo, M., Duarte, E., Marigliano, O., Görgen, C., Maestre, F. T. and Eisenhauer, N. 2021. Global projections of the soil microbiome in the Anthropocene. – Global Ecol. Biogeogr. 30: 987–999.

Guisan, A. and Rahbek, C. 2011. SESAM – a new framework integrating macroecological and species distribution models for predicting spatio-temporal patterns of species assemblages. – J. Biogeogr. 38: 1433–1444.

Guisan, A. and Thuiller, W. 2005. Predicting species distribution: offering more than simple habitat models. – Ecol. Lett. 8: 993–1009.

Hastie, T., Tibshirani, R. and Wainwright, M. 2015. Statistical learning with sparsity: the lasso and generalizations. – Chapman & Hall/CRC.

Henseler, J. 2021. Composite-based structural equation modeling: analyzing latent and emergent variables. – Guilford Press.

Hirzel, A. H. and Le Lay, G. 2008. Habitat suitability modelling and niche theory. – J. Appl. Ecol. 45: 1372–1381.

Hortal, J., de Bello, F., Diniz-Filho, J. A. F., Lewinsohn, T. M., Lobo, J. M. and Ladle, R. J. 2015. Seven shortfalls that beset large-scale knowledge of biodiversity. – Annu. Rev. Ecol. Evol. Syst. 46: 523–549.

Hutchinson, G. E. 1957. Concluding remarks. population studies: animal ecology and demography. – Cold Spring Harb. Symp. Quant. Biol. 22: 415–427.

Isaac, N. J. B., Jarzyna, M. A., Keil, P., Dambly, L. I., Boersch-Supan, P. H., Browning, E., Freeman, S. N., Golding, N., Guillera-Arroita, G., Henrys, P. A., Jarvis, S., Lahoz-Monfort, J., Pagel, J., Pescott, O. L., Schmucki, R., Simmonds, E. G. and O'Hara, R. B. 2020. Data integration for large-scale models of species distributions. – Trends Ecol. Evol. 35: 56–67.

Kass, J. M., Fukaya, K., Thuiller, W. and Mori, A. S. 2024. Biodiversity modeling advances will improve predictions of nature's contributions to people. – Trends Ecol. Evol. 39: 338–348.

Koller, D. and Freideman, N. 2009. Probabilistic graphical models: principles and techniques. – The MIT Press.

Larsen, P. E., Field, D. and Gilbert, J. A. 2012. Predicting bacterial community assemblages using an artificial neural network approach. – Nat. Methods 9: 621–625.

Leathwick, J. R. and Austin, M. P. 2001. Competitive interactions between tree species in New Zealand's old-growth indigenous forests. – Ecology 82: 2560.

Lefcheck, J. S. 2016. piecewiseSEM: piecewise structural equation modelling in r for ecology, evolution and systematics. – Methods Ecol. Evol. 7: 573–579.
Le Guillarme, N. and Thuiller, W. 2022. TaxoNERD: deep neural models for the recognition of taxonomic entities in the ecological and evolutionary literature. – Methods Ecol. Evol. 13: 625–641.

Lotka, A. J. 1920. Analytical note on certain rhythmic relations in organic systems. – Proc. Natl Acad. Sci. USA 6: 410–415.

MacKenzie, D. I., Nichols, J. D., Andrew Royle, J., Pollock, K. H., Bailey, L. L. and Hines, J. E. 2006. Occupancy estimation and modeling: inferring patterns and dynamics of species occurrence. – Academic Press.

MacLulich, D. A. 1936. Fluctuations in numbers of varying hares. – Science 83: 162–162.

Maiorano, L., Montemaggiori, A., Ficetola, G. F., O'Connor, L. and Thuiller, W. 2020. Tetra-EU 1.0: A species-level trophic metaweb of European tetrapods. – Global Ecol. Biogeogr. 29: 1452–1457.

Meynard, C. N., Leroy, B. and Kaplan, D. M. 2019. Testing methods in species distribution modelling using virtual species: what have we learnt and what are we missing? – Ecography 42: 2021–2036.

Montesinos-Navarro, A., Estrada, A., Font, X., Matias, M. G., Meireles, C., Mendoza, M., Honrado, J. P., Prasad, H. D., Vicente, J. R. and Early, R. 2018. Correction: community structure informs species geographic distributions. – PLoS One 13: e0200556.

Norberg, A. et al. 2019. A comprehensive evaluation of predictive performance of 33 species distribution models at species and community levels. – Ecol. Monogr. 89: e01370.

O'Connor, L. M. J., Pollock, L. J., Braga, J., Ficetola, G. F., Maiorano, L., Martinez-Almoyna, C., Montemaggiori, A., Ohlmann, M. and Thuiller, W. 2020. Unveiling the food webs of tetrapods across Europe through the prism of the Eltonian niche. – J. Biogeogr. 47: 181–192.

O'Hara, R. B. and Sillanpää, M. J. 2009. A review of Bayesian variable selection methods: what, how and which. – Bayesian Anal. 4: 85–117.

Ovaskainen, O., Tikhonov, G., Norberg, A., Guillaume Blanchet, F., Duan, L., Dunson, D., Roslin, T. and Abrego, N. 2017. How to make more out of community data? A conceptual framework and its implementation as models and software. – Ecol. Lett. 20: 561–576.

Paquette, A. and Hargreaves, A. L. 2021. Biotic interactions are more often important at species' warm versus cool range edges. – Ecol. Lett. 24: 2427–2438.

Pearl, J. 2009. Causal inference in statistics: an overview. – Stat. Surv. 3: 96–146.

Pichler, M., Boreux, V., Klein, A.-M., Schleuning, M. and Hartig, F. 2020. Machine learning algorithms to infer trait-matching and predict species interactions in ecological networks. – Methods Ecol. Evol. 11: 281–293.

Poggiato, G., Münkemüller, T., Bystrova, D., Arbel, J., Clark, J. S. and Thuiller, W. 2021. On the interpretations of joint modeling in community ecology. – Trends Ecol. Evol. 36: 391–401.

Poggiato, G., Andréoletti, J., Pollock, L. J. and Thuiller, W. 2024. Data from: Integrating food webs in species distribution models can improve ecological niche estimation and predictions. – Github Repository, https://github.com/giopogg/webSDM/tree/main/publication.

Pollock, L. J., O'Connor, L. M. J., Mokany, K., Rosauer, D. F., Talluto, L. and Thuiller, W. 2020. Protecting biodiversity (in all its complexity): new models and methods. – Trends Ecol. Evol. 35: 1119–1128.

Pollock, L. J., Tingley, R., Morris, W. K., Golding, N., O'Hara, R. B., Parris, K. M., Vesk, P. A. and McCarthy, M. A. 2014. Understanding co-occurrence by modelling species simultaneously with a joint species distribution model (JSDM). – Methods Ecol. Evol. 5: 397–406.

Potapov, A. M. 2022. Multifunctionality of belowground food webs: resource, size and spatial energy channels. – Biol. Rev. Camb. Philos. Soc. 97: 1691–1711.

Pottier, J., Dubuis, A., Pellissier, L., Maiorano, L., Rossier, L., Randin, C. F., Vittoz, P. and Guisan, A. 2013. The accuracy of plant assemblage prediction from species distribution models varies along environmental gradients. – Global Ecol. Biogeogr. 22: 52–63.

Ramazi, P., Kunegel-Lion, M., Greiner, R. and Lewis, M. A. 2021. Exploiting the full potential of Bayesian networks in predictive ecology. – Methods Ecol. Evol. 12: 135–149.

Read, J., Pfahringer, B., Holmes, G. and Frank, E. 2021. Classifier chains: a review and perspectives. – J. Artif. Intell. Res. 70: 683–718.

Renner, I. W., Elith, J., Baddeley, A., Fithian, W., Hastie, T., Phillips, S. J., Popovic, G. and Warton, D. I. 2015. Point process models for presence-only analysis. – Methods Ecol. Evol. 6: 366–379.

Rosa, I. M. D. et al. 2017. Multiscale scenarios for nature futures. – Nat. Ecol. Evol. 1: 1416–1419.

Shipley, B. 2000. Cause and correlation in biology: a User's guide to path analysis, structural equations and causal inference with R. – Cambridge Univ. Press.

Soberón, J. 2007. Grinnellian and Eltonian niches and geographic distributions of species. – Ecol. Lett. 10: 1115–1123.

Staniczenko, P. P. A., Sivasubramaniam, P., Suttle, K. B. and Pearson, R. G. 2017. Linking macroecology and community ecology: refining predictions of species distributions using biotic interaction networks. – Ecol. Lett. 20: 693–707.

Staniczenko, P. P. A., Suttle, K. B. and Pearson, R. G. 2018. Negative biotic interactions drive predictions of distributions for species from a grassland community. – Biol. Lett. 14: 20180426.

Strydom, T., Catchen, M. D., Banville, F., Caron, D., Dansereau, G., Desjardins-Proulx, P., Forero, N., Higino, G., Mercier, B., Gonzalez, A., Gravel, D., Pollock, L. and Poisot, T. 2021. A roadmap toward predicting species interaction networks (across space and time). – Philos. Trans. R. Soc. B 376: 20210063.

Strydom, T., Bouskila, S., Banville, F., Barros, C., Caron, D., Farrell, M. J., Fortin, M.-J., Hemming, V., Mercier, B., Pollock, L. J., Runghen, R., Dalla Riva, G. V. and Poisot, T. 2022. Food web reconstruction through phylogenetic transfer of low-rank network representation. – Methods Ecol. Evol. 13: 2838–2849.

Thuiller, W., Guéguen, M., Bison, M., Duparc, A., Garel, M., Loison, A., Renaud, J. and Poggiato, G. 2018. Combining point-process and landscape vegetation models to predict large herbivore distributions in space and time – a case study of *Rupicapra rupicapra*. – Divers. Distrib. 24: 352–362.

Thuiller, W., Münkemüller, T., Lavergne, S., Mouillot, D., Mouquet, N., Schiffers, K. and Gravel, D. 2013. A road map for integrating eco-evolutionary processes into biodiversity models. – Ecol. Lett. 16: 94–105.

Thuiller, W., Pollock, L. J., Gueguen, M. and Münkemüller, T. 2015. From species distributions to meta-communities. – Ecol. Lett. 18: 1321–1328.

Tredennick, A. T., Hooker, G., Ellner, S. P. and Adler, P. B. 2021. A practical guide to selecting models for exploration, inference and prediction in ecology. – Ecology 102: e03336.
van der Merwe, S., Greve, M., Olivier, B. and le Roux, P. C. 2021. Testing the role of functional trait expression in plant-plant facilitation. – Funct. Ecol. 35: 255–265.

Vehtari, A., Gelman, A. and Gabry, J. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. – Stat. Comput. 27: 1413–1432.

Volterra, V. 1926. Fluctuations in the abundance of a species considered mathematically. – Nature 118: 558–560.

Warton, D. I., Foster, S. D., De'ath, G., Stoklosa, J. and Dunstan, P. K. 2015a. Model-based thinking for community ecology. – Plant Ecol. 216: 669–682.

Warton, D. I., Blanchet, F. G., O'Hara, R. B., Ovaskainen, O., Taskinen, S., Walker, S. C. and Hui, F. K. C. 2015b. So many variables: joint modeling in community ecology. – Trends Ecol. Evol. 30: 766–779.

Windsor, F. M., van den Hoogen, J., Crowther, T. W. and Evans, D. M. 2022. Using ecological networks to answer questions in global biogeography and ecology. – J. Biogeogr. 50: 57–69.

Wisz, M. S. et al. 2013. The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling. – Biol. Rev. Camb. Philos. Soc. 88: 15–30.

Zurell, D., Berger, U., Cabral, J. S., Jeltsch, F., Meynard, C. N., Münkemüller, T., Nehrbass, N., Pagel, J., Reineking, B., Schröder, B. and Grimm, V. 2010. The virtual ecologist approach: simulating data and observers. – Oikos 119: 622–635.

Zurell, D., Pollock, L. J. and Thuiller, W. 2018. Do joint species distribution models reliably detect interspecific interactions from co-occurrence data in homogenous environments? – Ecography 41: 1812–1819.