# Modeling the effects of species associations and abiotic parameters on the abundance of mosquito species in a Mediterranean wetland

Ridwan Adeyemi Shittu, Stephanie Margarete Thomas, David A Roiz, Santiago Ruiz, Jordi Figuerola, Carl Beierkuhnlein

► To cite this version:

Ridwan Adeyemi Shittu, Stephanie Margarete Thomas, David A Roiz, Santiago Ruiz, Jordi Figuerola, et al.. Modeling the effects of species associations and abiotic parameters on the abundance of mosquito species in a Mediterranean wetland. Wetlands Ecology and Management, 2024, 32 (3), pp.381-395. 10.1007/s11273-024-09982-3 . hal-04806934

HAL Id: hal-04806934

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

Submitted on 27 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling the effects of species associations and abiotic parameters on the abundance of mosquito species in a Mediterranean wetland

Ridwan Adeyemi Shittu · Stephanie Margarete Thomas · David Roiz · Santiago Ruiz · Jordi Figuerola · Carl Beierkuhnlein

Received: 8 November 2023 / Accepted: 31 January 2024 / Published online: 11 March 2024
© The Author(s) 2024

**Abstract** Recent outbreaks of the West Nile virus have been reported in southern Spain, a region with important wetland habitats for migratory birds. Here, we analyzed the role of species association and abiotic parameters on the abundance of seven mosquito species in the Doñana National Park, Spain. We applied the Joint Species Distribution Models (JSDM), using the Hierarchical Modeling of Species Communities approach to simultaneously model the effect of habitat type, normalized difference vegetation index, hydroperiod, distance to rivers, land surface

temperature, and the abundance of seven mosquito species. We created three models with varied parameters and evaluated the effects of abiotic parameters and species-to-species associations, which served as proxies for species interactions. Our models estimated the species-specific responses of the mosquito community to abiotic parameters. A positive association was identified within the Culex community with potential biotic interactions. However, Ochlerotatus caspius, Ochlerotatus detritus, and Anopheles atroparvus had no statistically significant association with each other and with the Culex species after accounting for the effect of abiotic parameters. We also found that the potential species associations estimated at

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s11273-024-09982-3.

R. A. Shittu · S. M. Thomas (✉) · C. Beierkuhnlein
Department of Biogeography, University of Bayreuth,
Bayreuth, Germany
e-mail: stephanie.thomas@uni-bayreuth.de

R. A. Shittu
e-mail: ridwan.shittu@uni-bayreuth.de

C. Beierkuhnlein
e-mail: carl.beierkuhnlein@uni-bayreuth.de

S. M. Thomas · C. Beierkuhnlein
Center of Ecology and Environmental Research
(BayCEER), University of Bayreuth, Bayreuth, Germany

D. Roiz
MIVEGEC, IRD, CNRS, University of Montpellier,
Montpellier, France
e-mail: davidroiz@gmail.com

S. Ruiz
Servicio de Control de Mosquitos de La Diputación de
Huelva, Huelva, Spain
e-mail: sruiz@diphuelva.org

S. Ruiz · J. Figuerola
CIBER Epidemiology and Public Health - CIBERESP,
Madrid, Spain

J. Figuerola
Estación Biológica de Doñana - CSIC, Seville, Spain
e-mail: jordi@ebd.csic.es

C. Beierkuhnlein
Geographical Institute Bayreuth (GIB), University
of Bayreuth, Bayreuth, Germany
particular sites and abiotic parameters influenced the model projection for the average abundance of mosquito species. The JSDM will allow spatial projection of the abundance of each mosquito species, which is an important parameter for epidemiological models. The JSDM inference of species association is important because predation, competition, and facilitation affect the distribution and abundance of different species.

**Keywords** Biotic interactions · Species associations · Co-occurrence · JSDM · Mosquito community · Mosquito control · Species distribution model · Vector-borne diseases

## Background

Mosquito-borne diseases remain a global threat to human health and are re-emerging in Europe (Johnson et al. 2018). Climate change, biodiversity loss, global connectivity through travel and trade, and landscape changes, are considered the main drivers for expanding the distribution and abundance of mosquito species and the associated re-emergence of mosquito-borne diseases (Johnson et al. 2018; Bouzid et al. 2014; Tomasello and Schlagenhauf 2013).

Understanding the influence of abiotic and biotic conditions on the distribution and abundance of mosquitoes is a prerequisite for estimating the areas at risk of mosquito-borne diseases and establishing a reliable surveillance system (Golding et al. 2015). Many statistical approaches in this field only consider the effect of environmental variables and fail to consider the significant impact of biotic interactions among species, including facilitation, competition, and predation, and therefore overlook biodiversity effects. Utilizing Joint Species Distribution Models (JSDM) can reduce this gap and enhance research outcomes (Golding et al. 2015). Modeling the role of biotic interactions/associations, in addition to environmental niche modeling for individual species, enables quantifying this important contribution to mosquito species abundance and dynamics.

Biodiversity plays a complex role in disease transmission, particularly for vector-borne diseases (Levi et al. 2016). The diversity of species involved, including hosts, vectors, and pathogens, influences the dynamics and risks of these diseases. Two main

hypotheses, the “dilution effect” and the “amplification effect” may explain the relationship between biodiversity and disease transmission (Ostfeld and Keesing 2000; Keesing et al. 2006; Keesing and Ostfeld 2021). The dilution effect suggests that higher species diversity reduces pathogen transmission, resulting in a negative relationship with disease transmission. Conversely, the amplification effect describes situations where greater species diversity promotes pathogen transmission, leading to a positive relationship with disease risk (Keesing et al. 2006). The prevailing process depends on the spatial scale and context and is idiosyncratic (Salkeld et al. 2013). Changes in the abundance of host and vector species affect the likelihood of transmission of diseases like Lyme disease, West Nile fever, and the Usutu virus (Levi et al. 2016; Ostfeld and Keesing 2000; Johnson and Thieltges 2010; Roiz et al. 2019; Ferraguti et al. 2021).

The global and regional dynamics of the abundance and diversity of mosquitoes and mosquito-borne diseases are mainly due to changes in abiotic and biotic conditions (Burkett-Cadena et al. 2013; Ferraguti et al. 2016).

Burkett-Cadena et al. (Burkett-Cadena et al. 2013) linked mosquito abundance to host density, emphasizing the crucial role of host distribution in shaping vector patterns, reflecting parasitic symbiosis. Ferraguti et al. (Ferraguti et al. 2016) found higher mosquito abundance and species richness in natural and rural areas compared to urban settings, with links to environmental factors like land use, vegetation, and hydrology. Studies on the abundance of mosquito species have mostly focused on the implications of abiotic parameters with little focus on the potential influence of biotic parameters (Roiz et al. 2015; Ferraguti et al. 2021).

In (Roiz et al. 2015), the authors use remote sensing to analyze how environmental factors influence the distribution of seven mosquito species, including West Nile virus vectors. Roiz et al. (Roiz et al. 2015) findings highlight species-specific relationships, including a negative association between hydroperiod- as the number of days per year that an area of land is wet or the length of time that there is standing water at a location and mosquito presence, a positive correlation between Culex abundance and hydroperiod, and a positive link between NDVI and mosquito diversity. Inundation surfaces positively influence the mosquito species presence and mosquito richness,
excluding salt marsh mosquitoes. The study emphasizes landscape influence on mosquito distribution, pointing to potential changes in disease risk due to human landscape modification.

Most species outlined in Roiz et al.'s work (Roiz et al. 2015) have been documented as competent vectors of the West Nile virus (WNV). Those species identified as potential WNV vectors through a combination of vector competence experiments or field-based mosquito investigations include *Culex theileri*, *Culex modestus*, *Culex pipiens*, *Culex perexiguus*, *Ochlerotatus caspius*, and *Ochlerotatus detritus* (Jupp et al. 1972; McIntosh et al. 1976; Balenghien et al. 2007, 2008; Vázquez et al. 2011; Papa et al. 2011; Engler et al. 2013; Brustolin et al. 2016; Blagrove et al. 2016; Vogels et al. 2017; Mancini 2017; Kampen and Walther 2018).

In this study we applied the emerging approach of the Joint species distribution model (JSDM), specifically the Hierarchical Modeling of Species Communities (HMSC) to an existing dataset of the abundance of adult female mosquitoes (Roiz et al. 2015) to reanalyze the contribution of abiotic parameters. In addition, the potential species associations of mosquito species were estimated in a Mediterranean wetland. We estimated the taxa of mosquitoes that co-occur as a proxy of species interaction in mosquito communities. We expect that (a) abiotic parameters strongly influence the abundance of mosquitoes, (b) the random selection of the spatial location of traps on the different habitat types does not influence the abundance of each mosquito species, and (c) all mosquito species have potential biotic interaction after considering the impact of the abiotic parameters.

## Methods

### Study area and design

The study area was located at the Doñana Natural Space in Andalusia on the South-Western coast of Spain, covering an area of 106,000 ha. Doñana Natural space is a touristic designation but also serves as an important habitat for many species, including mosquitoes and migratory birds (Serrano et al. 2006; Green et al. 2016). About 5498.4 ha were sampled for mosquito species using BG-Sentinel traps from March to November 2010 (Roiz et al. 2015). The

sampling area was classified into six ecological units, and a stratified sampling design was applied, with the number of trapping sites in each ecological unit being proportional to its relative surface in the study area: marshland (22,109 ha, 47 traps), scrubland (11,582 ha, 31 traps), sand-dunes (9113 ha, 15 traps), fishpond (5473 ha, 10 traps), cropland (5207 ha, 6 traps) and rice fields (36,600 ha, 3 traps). There were fewer traps on the rice field because of theft and because the study was centered in the natural protected areas. We removed the traps on ricefields from analysis as the ricefields were subjected to control measures for mosquitoes during the sampling period. Though randomly located in different habitats, each trap serves as a sampling unit (Fig. 1).

### Mosquito abundance data

We extracted the relative abundance dataset of each mosquito species provided as a supplementary document in the publication by Roiz et al. (Roiz et al. 2015). The sampling units (trap locations) yielded mean abundance data for seven mosquito species—*Anopheles atroparvus*, *Culex perexiguus*, *Culex pipiens*, *Culex modestus*, *Culex theileri*, *Ochlerotatus caspius*, and *Ochlerotatus detritus*. Mosquito abundances was calculated as the mean of the number of females captured per trap night. A sampling unit that did not have abundance records for any of the species was removed. Therefore, a total of 108 sampling units were used for the analyses. The matrix of the mean annual abundance data of each of the seven most common mosquito species found in the wetland serves as the response variable (Ovaskainen et al. 2017).

### Abiotic parameters

The abiotic parameters describing the wetland were derived using remote sensing approaches. Remote sensing products were downloaded and analyzed (Landsat-5,7 image courtesy of the U.S. Geological Survey) from March to November 2010. The abiotic parameters include both continuous and categorical variables. The continuous abiotic parameters included the mean annual hydroperiod, the mean annual normalized difference vegetation index (NDVI), and the mean annual land surface temperature at trap sites. Parameters like NDVI and hydroperiod were derived
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Map of the Doñana wetlands. The white circle indicates the location of the traps. Adapted from Roiz et al., 2015</figcaption></figure>

using different buffer sizes (100 m, 250 m, 500 m, 1000 m, and 2000 m) around the 108 traps. We included the habitat types based on vegetation maps on which the traps were located as categorical parameters. Further methodological information on abiotic parameter calculations is described in Roiz et al. (Roiz et al. 2015). In addition, we calculated the land surface temperature (LST) for trap locations between March to November 2010 using the method described in Li et al. (2004) from remotely sensed data ( Landsat, https://earthexplorer.usgs.gov/).

## Model and model fitting

We used the variance inflation factor method to check for and reduce the effect of multicollinearity by removing the abiotic parameters that are highly correlated. Multicollinearity occurs when two or more

independent variables have a high correlation with one another in a regression model. The problem of multicollinearity makes it difficult to distinguish the effects of abiotic parameters from one another (Zuur et al. 2010; Cheng et al. 2022).

Joint species distribution models simultaneously estimate dependent variables’ responses to the independent variable(s) and evaluate species co-occurrence from residual covariance structure (Pollock et al. 2014). JSDM can concurrently estimate the effect of abiotic conditions on the abundance of each mosquito species and infer species-to-species associations from residual correlations as a proxy for biotic interactions (Golding et al. 2015; Pollock et al. 2014; Warton et al. 2015). JSDM explains variation in community assemblages by integrating contributions from the abiotic, biotic relationships, and perhaps spatially structured residual covariance (Pichler and
Hartig 2021). When using JSDMs, it is feasible to consider species distribution and community assemblage processes which is not possible with single-species distribution models (Wilkinson et al. 2021). Species assemblages can emerge from biotic interactions or responses to missing abiotic parameters (Golding et al. 2015; Pollock et al. 2014; Wilkinson et al. 2021).

We analysed the mosquito abundance data with Hierarchical Modelling of Species Communities (HMSC) (Ovaskainen et al. 2017; Ovaskainen and Abrego 2020). The HMSC model is a multivariate hierarchical generalized linear mixed model that includes both fixed and random effects. The unique extension of GLMM for JSDMs approaches like HMSC is the inclusion of latent variables, which can be described as unobserved predictors from missing parameters or biotic associations (Warton et al. 2015; Ovaskainen and Abrego 2020; Zhang et al. 2018).

The equations for the HMSC models employed in our analysis are outlined in Appendix S3. We used lognormal distribution due to the fact that the data of the mosquito abundance as the response variables (mean relative abundance) is count data with about 20% zero and not normal distribution (Ovaskainen and Abrego 2020). HMSC allows the modeling of mosquito species simultaneously as a response to abiotic conditions (fixed effect) while accounting for the spatial structure of the traps (random effect) in the Doñana natural space. Assuming the default prior distributions for the different model parameters detailed in Ovaskainen & Abrego, (Ovaskainen and Abrego 2020). We fitted the relative abundance of the mosquito species response to abiotic parameters and estimated potential species association using the HMSC model with the R-Package Hmsc (Tikhonov et al. 2020).

We fitted three models: the full model (FM), the environmental model (EM), and the spatial model (SM). The FM included both the selected abiotic parameters and the spatial random effect of the traps. The EM included only the selected abiotic parameters and no spatial random effect. The SM included only the spatial random effects of the traps and no abiotic parameters. We created the FM to estimate the influence of both the spatial random location of the traps (used to capture potential biotic interactions) and the selected abiotic parameters, while the EM and SM were created to assess the influence of the abiotic

parameters or the spatial random effect respectively on the abundance of each mosquito species. Species associations are estimated in the HMSC model’s random effect section which is structured using latent variables to capture the residual covariance estimated for the individual species. The potential biotic interactions/associations were modeled by including the random effect defined at the trap level. The trap-level random effect is utilized to pinpoint species associations that happen more or less frequently than would be predicted by chance. All the other HMSC parameters are the same for the models that are estimated.

Each HMSC model is fitted with Bayesian inference through the sampleMcmc function in the Hmsc package, and parameter inferences are estimated using the Markov Chain Monte Carlo (MCMC) approach because of the complex nature of the models (Golding et al. 2015). We tested different samples and thinning sizes to get the best for fitting the model. Thereafter, using the best-observed number of samples and thinning size, we sampled the posterior distribution for the model parameters using four MCMC chains, each of which was run for 1,000,000 iterations, of which the first 500,000 were removed as burn-in. The four chains were thinned by 1000 to yield 1000 posterior samples per chain and so 4000 posterior samples in total were used for further analysis. We examined the MCMC convergence by examining Gelman’s plot of the potential scale reduction factors (PSRF) of the model parameters. A PSRF value less than 1.2 is often recommended as representing satisfactory MCMC convergence as a larger PSRF (> 1.2) indicates a lack of convergence (Gelman and Rubin 1992).

## Model evaluation

We evaluated the models using explanatory and predictive power (Ovaskainen and Abrego 2020). The explanatory power, on the one hand, measures model performance by examining the accuracy of a model in predicting data used for model fitting. For the explanatory power, model predictions were based on the model fitted to all data. The predictive power, on the other hand, measures the effectiveness of a model to forecast data not used for model fitting. For this study, we calculated the predictive power using two-fold cross-validation, in which the sampling units were assigned randomly to two folds, and predictions
for each fold were based on the model fitted to data on the remaining fold. HMSC currently includes the cross-validation approach to evaluate the predictive power for model with lognormal function.

For the lognormal distribution model, both explanatory and predictive were measured by the coefficient of determination (R²) (Ovaskainen and Abrego 2020; Planillo et al. 2021). The coefficient of determination measures the proportion of variation of a dependent variable explained by an independent variable and remains an important metric for regression models (Wright 1921; Chicco et al. 2021). The closer to 1 for the values of R², the higher the proportion of the dependent variables that are explained by the independent variables.

We explored each mosquito species' response to the abiotic parameters considered in the FM and in the EM using variance partitioning analysis. The variance partitioning gives information on the importance of the abiotic parameters and potential biotic interaction in explaining the abundance of mosquito species in the wetland which help highlight the difference between the two models. Furthermore, we examined species-specific responses (species niches) estimated from the parameters of the models and assessed the pattern of the response whether positive, negative, or no significant response.

Similarly, we examined the species-to-species association revealed by the random effects included in the SM and FM. Raw association pertains to an association that can occur randomly when data is reorganized, whereas residual association signifies a form of environmentally conditioned analysis that factors in fluctuations linked to the alignment of species niches and environmental factors. Raw association corresponds to the outcome of SM, while Residual association corresponds to the outcome of FM. For the FM, positive species association suggests that species-to-species association occurs more frequently than would be expected by chance, whereas negative species association means the opposite (Tikhonov et al. 2020).

We used the fitted FM to make predictions that illustrate how the total mosquito species and each mosquito species abundance in the Doñana wetland are predicted to respond to changes in each of the abiotic parameters. We predicted the average total count of each mosquito species using the abiotic parameters used for the FM. For the prediction, two categories

of effects were estimated: total and marginal effects. The total effect measures the influence of one abiotic parameter (focus parameter) on the total count of all the mosquito species (or each mosquito species) when the remaining abiotic parameters (non-foci parameter) are at the most likely values for the values of the focus parameters. The most likely values of the other abiotic non-foci parameters were estimated using a generalized linear model (Ovaskainen and Abrego 2020). The marginal effect measures the influence of an abiotic parameter on the total count of all the mosquito species (or each mosquito species) when the remaining parameters are fixed at the mean values for the non-habitat and at the mode for habitat parameters. All analyses were performed in R version 4.1.2 (R Development Core Team, 2021).

## Results

We analyzed abundance data of the seven most common mosquito species found in the sampled area of the Doñana wetland. We found that *Oc. detritus* was absent in 40% of the sites. *Cx. theileri* was present in all the traps, and five other species (*Anopheles atroparvus*, *Oc. caspius*, *Cx. pipiens*, *Cx. modestus*, *Cx. perexiguus*) were captured in more than 70% of the sites. The descriptive information on the abundance of each mosquito species is reported in Appendix S1: Table S1.

The MCMC convergence of the HMSC models was found satisfactory using the potential scale reduction factors for the parameters of FM and EM. The PSRF was on average 1.02 for the FM and 1.00 for the EM. The Gelman's PSRF plot for the parameters of the FM are illustrated in Appendix S1: Figure S1.

The Full Model (FM) demonstrated a strong fit to the data, with mean R² values of 0.66 for explanatory power and 0.32 for predictive power. This indicates that the FM could account for approximately 66% of the variation observed in the abundance data for each mosquito species. Similarly, the Environmental Model (EM) exhibited a satisfactory model fit, with mean R² values of 0.46 for explanatory power and 0.32 for predictive power, signifying an ability to explain around 46% of the variation in the mosquito dataset.

Variance partitioning of both HMSC models with abiotic parameters as fixed effects explained a
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Variance partitioning among the abiotic parameters for mosquito abundance. The contribution of each variance component is given as a proportion of the total variance explained by the model given $R^2=46\%$, includes only abiotic parameters. **A** The environmental model (EM) average $R^2=66\%$, includes abiotic parameters (fixed effects) and potential biotic parameters captured at spatial locations of the traps (random effects)</figcaption></figure>

substantial amount of variation in the data of the average abundance of each mosquito species (Fig. 2). For the EM, distance to the river, annual hydroperiods, and habitat types explain about 25.8%, 25.8%, and 20.8% respectively of the variation in the abundance of mosquito data (Fig. 2A). The potential biotic interactions captured with the random effect of the spatial location of the traps contributed about 30.8% to the variation in the FM. (Fig. 2B). We found the distance to river, annual hydroperiod (buffered at 250 m) and habitat types explained about 18.9%, 15.9%, and 16.3% respectively, of the total variation explained by the FM (Fig. 2B). We found that the FM revealed the importance of the potential biotic interaction estimated in the abundance of each mosquito species.

We assessed the variation among the species by visualizing the species-specific response (species niches) using the parameters for two models that included the abiotic parameters (Fig. 3A and B). The abiotic parameters, affecting each mosquito species differ significantly with at least 0.95 posterior probability.

For the EM, the habitat preferences of the species were in comparison to cropland which serves as the intercept by the HMSC model. Most of the studied mosquito species showed positive preference for cropland except *Oc. detritus* and *An. atroparvus*. This model also showed statistical support for the positive preference of sand dune for the abundance of *An. atroparvus*, *Cx. modestus*, *Cx. perexiguus*, *Cx*,

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Heatmap of estimated species niches: species responses to abiotic parameters. **A** The environmental model (EM). **B** The full model (FM). The red and blue colors show the parameters estimated to be positive and negative responses, respectively, with at least 95% posterior probability. The white color</figcaption></figure>

shows responses that did not have statistical support. The estimated niches captured by both models (EM and FM) are depicted in darker colors, while the lighter colors illustrate the estimated niches unique to the EM model
pipiens, and *Cx. theileri*. *Cx. theileri*, *Cx. modestus* and *Oc. detritus* showed a positive preference for fishponds. *Cx. theileri*, and *Oc. caspius* showed positive preference for marshlands, while *Oc. detritus* showed a negative preference. Similarly, *Oc. caspius* and *Cx. theileri* showed positive preference for scrublands.

The continuous variable of the EM showed different species-specific responses (Fig. 3A). The distance to river was related to the annual abundance of each of the mosquito species except for *Cx. perexiguus*. The abundance of *Cx. theileri*, *Cx. pipiens* and *Cx. modestus* showed a positive relationship, while *An. atroparvus*, *Oc. caspius* and *Oc. detritus* showed a negative relationship, with annual hydroperiod. Annual NDVI was positively related to the abundance of *Cx. perexiguus*, *Cx. modestus* and *Cx. pipiens*. The distance to rice fields was positively related to *Cx. pipiens* and *Cx. perexiguus* respectively. Only *Cx. theileri* abundance was positively related to distance to urban areas. The annual mean land surface temperature was positively related to the abundance of *An. atroparvus*, *Cx. perexiguus* and *Cx. pipiens*.

We found changes in the species-specific response for the FM which controls for the effect of the random spatial location of the traps compared to the EM. For the habitats, *Cx. theileri* and *Cx. perexiguus* showed a positive preference for cropland, the other species shows no statistically significant preference for cropland in comparison to the EM. Many of the positive and negative responses of the mosquito species in the EM were not statistically significant after controlling for the spatial random effect of traps (Fig. 3B). However, hydroperiod, distance to river and habitats were found by both EM and FM to significantly influence the abundance of the mosquito species in the wetland.

The preferences of the species are further explored using the prediction output of the fitted FM. We found that the mosquitoes abundance are predicted to vary to changing abiotic parameters. The FM predicted that Sand dunes have a positive relationship with the total mosquito abundance than the other habitats (Fig. 4A). However, *Cx. theileri*, the most abundant species was found and predicted to be more abundant in marshland than in Sand dunes (Fig. 4B), while the *Oc. detritus* being the least abundant species was found and predicted to be more abundant in fishponds than sand dunes (Fig. 4C). The estimated model prediction of the total and marginal effect of abiotic parameters on the remaining mosquito species

over abiotic parameters (categorical and continuous) are illustrated in Appendix S2.

The SM shows raw associations of the species without accounting for the effect of abiotic parameters while the FM shows residual associations for the species after accounting for the effect of abiotic parameters. We found different associations for each mosquito species in both FM and SM (model that included random effects of traps in the model formulation). The SM estimated a positive raw association between *Cx. theileri*, *Cx. pipiens*, and *Cx. modestus*, a negative raw association between *Cx. modestus* and with *Oc. detritus* (Fig. 5A). Likewise, a positive raw association was found among *Cx. pipiens*, *Cx. perexiguus* and *An. atroparvus*. The association of species for the FM differs from the SM. The FM estimated a positive pattern of associations for the *Culex* species; the other species showed no significant association after controlling for the effect of abiotic parameters (Fig. 5B).

## Discussion

Understanding the factors influencing the distribution of mosquito communities is important to enhance the reduction in uncertainties and challenges in mosquito-borne disease risks and transmission. Our analysis supports that both abiotic and biotic factors influenced the distribution and abundance of the seven most common mosquito species in the Doñana wetland. The potential biotic associations are estimated to be positive for the *Culex* species found in Doñana.

Both the present study and the study conducted by Roiz et al. (2015) share several similarities in their examination of the relationship between abiotic parameters and mosquito abundance, analysing annual relative abundance by utilizing data from the same source (Roiz et al. 2015). Both studies identify hydroperiod as a significant factor influencing mosquito abundance. Specifically, they find a positive relationship between hydroperiod and the abundance of *Cx. theileri*, *Cx. pipiens*, and *Cx. modestus*. The studies also found a positive association between NDVI and the annual abundance of *Cx. perexiguus* and *Cx. pipiens*.

However, there are notable differences between the two studies in their approaches and some of their findings. The modeling approaches employed differ, with
**Fig. 4** Predictions of species abundance over the habitat types in the full model (FM). The predictions have been made by varying the landscapes and setting the other abiotic parameters to their expectation condition on each habitat. **A** Total mosquito abundance, **B** *Cx. theileri* abundance, and **C** *Oc. detritus* abundance. The figures are based on the total effect estimation. The ellipses on the graph represent the distribution of the abundance data of mosquito. The error bar represents the credible interval of the posterior distribution (similar to standard error of mean).

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

<figure><img src="image_6.png" /></figure>
<caption>**Fig. 4** (continued</caption>

C

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Residual species to species associations for mosquito abundance. **A** The space model (SM), **B** The full model (FM). The red and blue color indicates those species pairs with at least 95% posterior probability support for positive and negative associations respectively</figcaption></figure>

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

Roiz et al. (2015) using generalized linear models estimating model parameters using maximum likelihood approach, while the present study utilizes an extended generalized linear models approach termed JSDM built on a Bayesian framework. Additionally, the specific parameters used in the models also differ between the two studies. The present study found a positive relationship between hydroperiod and the

annual abundance of *Cx. perexiguus*. This study identifies negative associations between hydroperiod and the annual abundance of *Oc. detritus*, *Oc. caspius*, and *An. atroparvus*.

We demonstrated that including other parameters like land surface temperature and estimation of potential biotic species associations using JSDM could help reveal important species association patterns
and improved the estimation of the influence of abiotic parameters for the mosquito species (EM-average R²=46% and FM-average R²=66%).

Similar to other previous modeling options, the models used in this study were constructed based on many assumptions (i.e., especially given the Bayesian framework that relies mainly on prior distribution for the data) which may not agree with the community ecology of the studied species (Ovaskainen and Abrego 2020). Some studies using similar JSDM have identified important ecological relationships between species. For example, Lany et al. (2018) analyzed the abundance of two invasive insects and concluded that incorporating spatial and temporal dependence using JSDM revealed the dependence of a species' abundance on other species in the community. Abrego et al. (2021) reported increased predictive power for the arctic arthropod species response to increasing temperature using JSDM, which included the potential species interaction rather than the model based solely on abiotic parameters. Nevertheless, it is important to check for the limitations of modeling approaches like JSDM and that proper findings might arise from combining model results with data from experiments and functional markers that are appropriate to the studied biotic interaction (Dormann et al. 2018).

We found statistical support for the influence of habitats and other abiotic parameters on the abundance of mosquito species in the wetland. The FM revealed that abiotic parameters with greater impact include habitat type (16.3%), hydroperiod (15.9%), and distance to the river (18.9%) compared to parameters such as land-surface temperature (4.4%), vegetation index (4.7%), and distance to ricefield (5.8%) or urban areas (3.2%).

Although the abundance of the mosquito species varied from one habitat to another, the abundance was higher in areas that might have a lower impact from humans such as the marshland, sand dunes, and scrubland than areas like the cropland except for fishponds. This observation agrees partially with a previous study (Ferraguti et al. 2016). Ferraguti et al. (2016) evaluated the effect of urbanization on mosquito communities and found that mosquitoes were more abundant in natural and rural areas than the urban regions.

We discovered that the potential biotic associations around the spatial distribution of the included

traps had an impact on the fluctuation in mosquito abundance within the wetland. Our full model estimated that about 30.8% of the variation in abundance is attributed to potential species-to-species associations.

The potential biotic interaction (captured in the form of species-to-species association) among the *Culex* species was higher than what would be expected by chance after controlling for the effect of abiotic parameters. Our analysis suggests that the potential biotic interaction estimated from JSDM might not only be on evolutionary similarity in the habitat requirements of the mosquito species in the wetland but also reflect a joint response to predation or other biotic or abiotic parameters not included in the model. The distribution patterns of mosquitoes are known to be influenced by various factors related to their habitats, including both natural characteristics and human activities. These factors encompass the physical and chemical properties of the habitats, such as pH and salinity, as well as anthropogenic influences like the presence of sewage and fertilizers (Multini et al. 2021; Neogi et al. 2014), that may habitat suitability of egg laying by female mosquitoes and/or larvae development. Notably, experimental studies have provided insights into the impact of salinity on the colonization behavior and composition of mosquito assemblages, particularly those species that are sensitive to changes in salinity levels (Neogi et al. 2014). While it is worth mentioning that these experimental findings may not align entirely with observations from field studies, they contribute valuable information toward understanding the role of salinity in shaping mosquito distribution patterns (Carver et al. 2009).

The study by Golding et al. (Golding et al. 2015) showed that ditch shrimp and fish are negatively linked with the larvae of *An. maculipennis* and *Cx. modestus* after controlling for abiotic parameters.

In addition, Ohba et al. (Ohba et al. 2012), reported that *Culex tritaeniorhynchus* females are more sensitive when choosing oviposition sites to predators' presence cues than *Aedes albopictus*. Consequently, *Cx. tritaeniorhynchus* utilize predator cues to choose the safe places to lay their eggs. As a result, traits like predator sensitivity may potentially affect other *Culex* species and results in the positive species association estimated by our model for the *Culex* species in the mosquito abundance dataset we analyzed.
The association between mosquito species indicated by our full-model analysis, may be influenced by the similarity in their feeding preferences. However, *Cx. theileri* is a mammophilic species while *Cx. perexiguus*, *Cx. modestus* and *Cx. pipiens* are mainly ornithophilic species (Muñoz et al. 2012; Gómez-Díaz and Figuerola 2010). For WNV vectors such as *Cx. perexiguus* and *Cx. pipiens*, they feed mainly in birds with strong preferences for some avian species, such as the European blackbirds, that is bitten at a higher rate than expected based on their abundance relative to other avian species (Rizzoli et al. 2015a). The presence of competent mosquito vectors, their feeding preferences, and the composition of the avian host community influences the transmission of pathogens (Ferraguti et al. 2021; Marm Kilpatrick et al. 2006).

Understanding the co-occurrence of mosquito species is essential for establishing an effective control system. Figuerola (2022) suggested that the proliferation of *Cx. perexiguus* in certain southern regions of Spain, facilitated by the absence of mosquito control measures, may have triggered the outbreak of WNV in the summer and spring of 2020 (Figuerola et al. 2022; Rodríguez-Alarcón et al. 2021). *Culex perexiguus* mosquitoes likely played a central role in the enzootic transmission of WNV in the surrounding areas and within some villages, while *Culex pipiens* may have acted as a bridge vector transmitting the virus to humans (Figuerola et al. 2022; Rizzoli et al. 2015b). In fact, Figuerola et al. (2022) suggested that human infections were more likely caused by mosquito bites within the villages rather than in natural areas or nearby rice fields. It is important to note that, so far, no human cases of WNV have been reported in the villages closer to the Doñana Natural Space, which is the study area.

The four species of *Culex* estimated as having potential biotic interaction with one another are all potential vectors of mosquito-borne diseases. Our analysis shows that potential positive biotic interactions among *Cx. species* might indicate that the abundance of *Cx. theileri* favours the abundance of *Cx. pipiens*. A conceivable interpretation for the observed intraspecific interaction is that the existence of larvae from a more prevalent species, such as *Culex theileri*, may mitigate predation on the less abundant species like *Culex pipiens*. Baldacchino et al. (Baldacchino et al. 2017) studied larval competition between

invasive mosquitoes *Aedes koreicus* and *Aedes albopictus* in northeastern Italy, observed weak competition between the two species and concluded that *Aedes albopictus* develops faster and positively influences larger *Ae. koreicus* emergence.

The positive correlation among the four *Culex* species could also be attributed to a collective response to unmeasured environmental variables or resources that remain unaccounted for in the Joint Species Distribution Model (JSDM) (Golding et al. 2015; Pollock et al. 2014). The co-occurrence of *Cx. pipiens* and *Cx. perexiguus* may favour West Nile virus amplification and transmission to humans because *Cx. perexiguus* is a key species for WNV amplification in southwestern Spain (Muñoz et al. 2012; Figuerola et al. 2022; Ferraguti et al. 2021).

Modeling a few of the species found in a similar habitat together with mosquitoes could reveal important species that are associated with mosquitoes. The pattern of the association that could be deduced through modeling might help enhance the possibilities of using biological methods for controlling mosquito distribution and abundance. Therefore, expanding the possibility of reducing the burden of mosquito-borne diseases like WNV. In addition, results from estimating biotic interaction could be used in epidemiological models, and hybrid models (that combine features from both mechanistic and correlative modeling approaches) to fully understand the structure and dynamics of relationships among mosquitoes, hosts, pathogens, and their habitats.

## Conclusion

This study concludes that potential species associations and abiotic parameters influence the abundance of mosquito communities in the Doñana wetlands using the joint species distribution modeling approach. The influence of abiotic parameters was more than the potential biotic interactions. Estimated species-to-species association among *Culex* species revealed the potential of interactions between the species, which are the vector of mosquito-borne diseases like West Nile virus, commonly found in the southern part of Spain. We suggest including other species (predators, host species via blood meals, plant species male mosquito feed on) commonly found in similar habitats to mosquitoes in joint species distribution
modeling studies could further enhance the detection of biotic interactions among mosquitoes and other flora and fauna. Estimating interactions/associations among different species in similar habitats with mosquitoes could enhance the understanding of the risk of emergence of mosquito-borne pathogens and the influence of biodiversity, with the aim to develop proactive control actions.

**Acknowledgements** We want to thank Osto Ovaskainen's group for organizing a workshop that facilitate the use of JSDM specifically the HMSC approach. Maximilian Pichler for his assistance in the research model development. Their expertise and guidance significantly enriched our study, and we are truly appreciative of their support.

**Author contributions** SMT and RAS conceptualized the study, wrote the original draft, and reviewed and edited the manuscript. DR, SR, and JF provided the dataset. RS performed the data analysis and visualized the results. SMT and CB supervised the study and acquired the funding. SMT was responsible for project management. All authors reviewed and edited the manuscript and have read and agreed to the published version of the manuscript.

**Funding** Open Access funding enabled and organized by Projekt DEAL. This research was funded through the 2018–2019 BiodivERsA joint call for research proposals, under the BiodivERsA3 ERA-Net COFUND program, and with the funding organization DFG, German Research Foundation (BE 2192/18-1). This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—491183248.

**Data availability** The datasets analyzed for this study can be obtained from Roiz et al. 2015 or by contacting the corresponding author.

## Declarations

**Competing interests** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly

from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Abrego N, Roslin T, Huotari T, Ji Y, Schmidt NM, Wang J et al (2021) Accounting for species interactions is necessary for predicting how arctic arthropod communities respond to climate change. *Ecography* 44:885–896

Baldacchino F, Arnoldi D, Lapère C, Rosà R, Montarsi F, Capelli G et al (2017) Weak larval competition between two invasive mosquitoes *Aedes koreicus* and *Aedes albopictus* (Diptera: Culicidae). *J Med Entomol* 54:1266–1272

Balenghien T, Vazeille M, Reiter P, Schaffner F, Zeller H, Bicout DJ (2007) Evidence of laboratory vector competence of culex modestus for west Nile virus. *J Am Mosq Control Assoc* 23:233–236

Balenghien T, Vazeille M, Grandadam M, Schaffner F, Zeller H, Reiter P et al (2008) Vector competence of some french culex and aedes mosquitoes for west Nile virus. *Vector-Borne Zoonotic Dis* 8:589–596

Blagrove MSC, Sherlock K, Chapman GE, Impoinvil DE, McCall PJ, Medlock JM et al (2016) Evaluation of the vector competence of a native UK mosquito Ochlerotatus detritus (*Aedes detritus*) for dengue, chikungunya and West Nile viruses. *Parasit Vectors* 9:452

Bouzid M, Colón-González FJ, Lung T, Lake IR, Hunter PR (2014) Climate change and the emergence of vector-borne diseases in Europe: case study of dengue fever. *BMC Public Health* 14:781

Brustolin M, Talavera S, Santamaría C, Rivas R, Pujol N, Aranda C et al (2016) *Culex pipiens* and *Stegomyia albopicta* (= *Aedes albopictus*) populations as vectors for lineage 1 and 2 West Nile virus in Europe. *Med Vet Entomol* 30:166–173

Burkett-Cadena ND, McClure CJW, Estep LK, Eubanks MD (2013) Hosts or habitats: what drives the spatial distribution of mosquitoes? *Ecosphere*. 4:art30

Carver S, Storey A, Spafford H, Lynas J, Chandler L, Weinstein P (2009) Salinity as a driver of aquatic invertebrate colonisation behaviour and distribution in the wheatbelt of Western Australia. *Hydrobiologia* 617:75–90

Cheng J, Sun J, Yao K, Xu M, Cao Y (2022) A variable selection method based on mutual information and variance inflation factor. *Spectrochim Acta A Mol Biomol Spectrosc* 268:120652

Chicco D, Warrens MJ, Jurman G (2021) The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation. *PeerJ Comput Sci* 7:e623

Dormann CF, Bobrowski M, Dehling DM, Harris DJ, Hartig F, Lischke H et al (2018) Biotic interactions in species distribution modelling: 10 questions to guide interpretation and avoid false conclusions. *Glob Ecol Biogeogr* 27:1004–1016

Engler O, Savini G, Papa A, Figuerola J, Groschup M, Kampen H et al (2013) European surveillance for west Nile virus
in mosquito populations. Int J Environ Res Public Health 10:4869–4895

Ferraguti M, La Martínez-de PJ, Roiz D, Ruiz S, Soriguer R, Figuerola J (2016) Effects of landscape anthropization on mosquito community composition and abundance. Sci Rep 6:29002

Ferraguti M, Heesterbeek H, Martínez-de La Puente J, Jiménez-Clavero MÁ, Vázquez A, Ruiz S et al (2021) The role of different *Culex* mosquito species in the transmission of West Nile virus and avian malaria parasites in Mediterranean areas. Transbound Emerg Dis 68:920–930

Ferraguti M, Martínez-de La Puente J, Figuerola J (2021) Ecological effects on the dynamics of west nile virus and avian plasmodium: the importance of mosquito communities and landscape. Viruses. 13:1208

Ferraguti M, Martínez-de La Puente J, Jiménez-Clavero MÁ, Llorente F, Roiz D, Ruiz S et al (2021) A field test of the dilution effect hypothesis in four avian multi-host pathogens Garcia-Arenal F., editor. PLoS Pathog. 17:e1009637

Figuerola J, Jiménez-Clavero MÁ, Ruíz-López MJ, Llorente F, Ruiz S, Hoefer A et al (2022) A one health view of the West Nile virus outbreak in Andalusia (Spain) in 2020. Emerg Microbes Infect 11:2570–2578

García San Miguel Rodríguez-Alarcón L, Fernández-Martínez B, Sierra Moros MJ, Vázquez A, Julián Pachés P, García Villacieros E et al (2021) Unprecedented increase of West Nile virus neuroinvasive disease Spain summer. Eurosurveillance. https://doi.org/10.2807/1560-7917.ES.2021.26.19.2002010

Gelman A, Rubin DB (1992) Inference from iterative simulation using multiple sequences. Stat Sci. https://doi.org/10.1214/ss/1177011136.full

Golding N, Nunn MA, Purse BV (2015) Identifying biotic interactions which drive the spatial distribution of a mosquito community. Parasit Vectors 8:367

Gómez-Díaz E, Figuerola J (2010) New perspectives in tracing vector-borne interaction networks. Trends Parasitol 26:470–476

Green AJ, Bustamante J, Janss GFE, Fernández-Zamudio R, Díaz-Paniagua C (2016) Doñana Wetlands (Spain). In: Finlayson CM, Milton GR, Prentice RC, Davidson NC (eds) Wetl Book. Springer, Dordrecht, pp 1–14

Johnson PTJ, Thieltgess DW (2010) Diversity, decoys and the dilution effect: how ecological communities affect disease risk. J Exp Biol 213:961–970

Johnson N, Fernández de Marco M, Giovannini A, Ippoliti C, Danzetta ML, Svartz G et al (2018) Emerging mosquito-borne threats and the response from European and eastern Mediterranean countries. Int J Environ Res Public Health. 15:2775

Jupp PG, McIntosh BM, Dickinson DB (1972) Quantitative experiments on the vector capability of culex (Culex) Theileri theobald with west nile and sindbis viruses1. J Med Entomol 9:393–395

Kampen H, Walther D (2018) Vector Potential of Mosquito Species (Diptera: Culicidae) Occurring in Central Europe. In: Dis Mosq-Borne (ed) Benelli G, Mehlhorn H. Cham, Springer, Cham, pp 41–68. https://doi.org/10.1007/978-3-319-94075-5_5

Keesing F, Ostfeld RS (2021) Dilution effects in disease ecology. Ecol Lett 24:2490–2505

Keesing F, Holt RD, Ostfeld RS (2006) Effects of species diversity on disease risk. Ecol Lett 9:485–498

Lany NK, Zarnetske PL, Schliep EM, Schaeffer RN, Orians CM, Orwig DA et al (2018) Asymmetric biotic interactions and abiotic niche differences revealed by a dynamic joint species distribution model. Ecology 99:1018–1023

Levi T, Keesing F, Holt RD, Barfield M, Ostfeld RS (2016) Quantifying dilution and amplification in a community of hosts for tick-borne pathogens. Ecol Appl 26:484–498

Li F (2004) Deriving land surface temperature from Landsat 5 and 7 during SMEX02/SMACEX. Remote Sens Environ 92:521–534

Mancini, (2017) Mosquito species involved in the circulation of West Nile and Usutu viruses in Italy. Vet Ital 53:97–110

Marm Kilpatrick A, Daszak P, Jones MJ, Marra PP, Kramer LD (2006) Host heterogeneity dominates West Nile virus transmission. Proc R Soc B 273:2327–2333

Mcintosh BM, Jupp PG, Dos-Santos I, Meenehan GM (1976) Epidemics of west nile and sindbis viruses in south Africa with Culex univittatus as vector. South Afr J Sci 72:295–300

Multini LC, Oliveira-Christe R, Medeiros-Sousa AR, Evangelista E, Barrio-Nuevo KM, Mucci LF et al (2021) The influence of the pH and salinity of water in breeding sites on the occurrence and community composition of immature mosquitoes in the green belt of the city of São Paulo. Brazil Insects 12:797

Muñoz J, Ruiz S, Soriguer R, Alcaide M, Viana DS, Roiz D et al (2012) Feeding patterns of potential west nile virus vectors in south-west Spain schneider BS, editor. PLoS ONE 7:e39549

Neogi SB, Yamasaki S, Alam M, Lara RJ (2014) The role of wetland microinvertebrates in spreading human diseases. Wetl Ecol Manag 22:469–491

Ohba S-Y, Ohtsuka M, Sunahara T, Sonoda Y, Kawashima E, Takagi M (2012) Differential responses to predator cues between two mosquito species breeding in different habitats. Ecol Entomol 37:410–418

Ostfeld RS, Keesing F (2000) Biodiversity and disease risk: the case of lyme disease. Conserv Biol 14:722–728

Ovaskainen O, Abrego N (2020) Joint species distribution modelling. Cambridge University Press, Cambridge

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

Papa A, Xanthopoulou K, Gewehr S, Mourelatos S (2011) Detection of West Nile virus lineage 2 in mosquitoes during a human outbreak in Greece. Clin Microbiol Infect 17:1176–1180

Pichler M, Hartig F (2021) A new joint species distribution model for faster and more accurate inference of species associations from big community data. Methods Ecol Evol 12:2159–2173

Planillo A, Kramer-Schadt S, Buchholz S, Gras P, Von Der Lippe M, Radchuk V (2021) Arthropod abundance modulates bird community responses to urbanization. Divers Distrib 27:34–49

Pollock LJ, Tingley R, Morris WK, Golding N, O’Hara RB, Parris KM et al (2014) Understanding co-occurrence by modelling species simultaneously with a joint species
distribution model (JSDM). McPherson J, editor. Methods Ecol Evol 5:397–406

Rizzoli A, Bolzoni L, Chadwick EA, Capelli G, Montarsi F, Grisenti M et al (2015a) Understanding West Nile virus ecology in Europe: *Culex pipiens* host feeding preference in a hotspot of virus emergence. *Parasit Vectors* 8:213

Rizzoli A, Jiménez-Clavero MA, Barzon L, Cordioli P, Figuerola J, Koraka P et al (2015b) The challenge of West Nile virus in Europe: knowledge gaps and research priorities. *Eurosurveillance*. https://doi.org/10.2807/1560-7917. ES2015.20.20.21135

Roiz D, Ruiz S, Soriguer R, Figuerola J (2015) Landscape effects on the presence, abundance and diversity of mosquitoes in Mediterranean wetlands. PLoS ONE 10:e0128112

Roiz D, Vázquez A, Ruiz S, Tenorio A, Soriguer R, Figuerola J (2019) Evidence that passerine birds act as amplifying hosts for usutu virus circulation. *EcoHealth* 16:734–742

Salkeld DJ, Padgett KA, Jones JH (2013) A meta-analysis suggesting that the relationship between biodiversity and risk of zoonotic pathogen transmission is idiosyncratic. *Laferty* K, editor. *Ecol Lett* 16:679–686

Serrano L, Reina M, Martín G, Reyes I, Arechederra A, León D et al (2006) The aquatic systems of Doñana (SW Spain): watersheds and frontiers. *Limnetica* 25:11–32

Tikhonov G, Opedal ØH, Abrego N, Lehikoinen A, de Jonge MMJ, Oksanen J et al (2020) Joint species distribution modelling with the r-package Hmsc. *Methods Ecol Evol* 11:442–447

Tomasello D, Schlagenhauf P (2013) Chikungunya and dengue autochthonous cases in Europe, 2007–2012. *Travel Med Infect Dis* 11:274–284

Vázquez A, Ruiz S, Herrero L, Moreno J, Molero F, Magallanes A et al (2011) West Nile and Usutu viruses in mosquitoes in Spain, 2008–2009. *Am J Trop Med Hyg* 85:178–181

Vogels CB, Göertz GP, Pijlman GP, Koenraadt CJ (2017) Vector competence of European mosquitoes for West Nile virus. *Emerg Microbes Infect* 6:1–13

Warton DI, Blanchet FG, O’Hara RB, Ovaskainen O, Taskinen S, Walker SC et al (2015) So Many variables: joint modeling in community ecology. *Trends Ecol Evol* 30:766–779

Wilkinson DP, Golding N, Guillera-Arroita G, Tingley R, McCarthy MA (2021) Defining and evaluating predictions of joint species distribution models. *Methods Ecol Evol* 12:394–404

Wright S (1921) Correlation and Causation. *J Agric Res* 20:557–585

Zhang C, Chen Y, Xu B, Xue Y, Ren Y (2018) Comparing the prediction of joint species distribution models with respect to characteristics of sampling data. *Ecography* 41:1876–1887

Zuur AF, Ieno EN, Elphick CS (2010) A protocol for data exploration to avoid common statistical problems: *Data exploration*. *Methods Ecol Evol* 1:3–14

**Publisher’s Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.