Refuge distributions and landscape connectivity affect host-parasitoid dynamics: Motivations for biological control in agroecosystems

Abstract

Species distributions are affected by landscape structure at different spatial scales. Here we study how the interplay between dispersal at different spatial scales and landscape connectivity and composition affect local species dynamics. Using a host-parasitoid model, we assessed host density and host occupancy on the landscape, under different parasitoid dispersal ranges and three local distributions of non-crop habitats, areas where hosts are unable to grow but parasitoids are provided with alternative hosts and food resources. Our results show distinct responses of host density to increases in non-crop area, measured by differences in slopes for different distributions of non-crop habitats, and that the effect of local landscape composition on species dynamics depends on the landscape connectivity at the regional scale. Moreover, we show how host density and occupancy are affected by increasing parasitoid dispersal ranges depending on landscape structure. Our results demonstrate the role of local and regional scales on species distributions and the importance of the combined effects of species biological parameters and landscape structure on species dynamics. Finally, we highlight the relevance of these aspects for the development of better strategies of biological control.

About the Authors:

Lucas D. Fernandes

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Department of Life Sciences, Imperial College London, Silwood Park, Ascot, Berkshire, United Kingdom, Departamento de Entomologia e Acarologia, Escola Superior de Agricultura Luiz de Queiroz—Universidade de São Paulo (USP), Piracicaba, SP, Brazil

ORCID logo https://orcid.org/0000-0001-6053-2277

Angelica S. Mata

Roles Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Writing – original draft, Writing – review & editing

Affiliation: Departamento de Física, Universidade Federal de Lavras (UFLA), Lavras, MG, Brazil

Wesley A. C. Godoy

Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing

Affiliation: Departamento de Entomologia e Acarologia, Escola Superior de Agricultura Luiz de Queiroz—Universidade de São Paulo (USP), Piracicaba, SP, Brazil

Carolina Reigada

Roles Conceptualization, Formal analysis, Methodology, Software, Supervision, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Departamento de Ecologia e Biologia Evolutiva, Universidade Federal de São Carlos (UFSCar), Rodovia Washington Luís, São Carlos, SP, Brazil

Introduction

The importance of landscape structure in influencing species distributions and metapopulation dynamics has been recognized by both theoretical and experimental studies [1–4]. Although it has been shown that heterogeneous landscape structure and specific dispersal patterns can alter population viability and species coexistence [5, 6], studies usually focus at a particular scale of spatial organization. Since there is now growing empirical evidence that species dynamics are in many cases subject to multi-scale processes [7–9], with important drivers from local to regional scales, the question of how these different scales simultaneously affect species distributions still needs to be resolved.

This question is particularly relevant for agriculture, as landscape composition and structure, by conditioning movement patterns, directly influence density and occupation of insects and potential agricultural pests [10–13]. Furthermore, the definition of biological control strategies strongly depends on how landscape connectivity and local resource quality influence the dynamics of interaction between pests and control agents [14, 15].

Hymenopterous parasitoids are important control agents of economically relevant pests in both agricultural and natural environments. These insects have foraging behaviours that vary according to landscape composition, which influences their distribution and consequently the regulation of host populations [16, 17]. Generalist parasitoids can also play an important role on pest management, as their presence in landscapes combining crop and non-crop areas may help controlling pest populations [18]. When the insect pest is not present or present in low densities on the crop, the existence of alternative plants can promote refuge habitat for parasitoids, allowing their persistence in the area and collaborating to the prevention of future pest irruptions [19].

Many studies have shown the positive effects of habitat heterogeneity on the increase of abundance and diversity of natural enemies in agricultural ecosystems [12, 20, 21]. In this context, the management of crop areas to increase the heterogeneity of agricultural landscapes is a viable alternative in conservation biological control programs. Non-crop habitats support predators and parasitoids providing alternative prey/hosts and food resources, shelter, and overwintering areas, while also facilitating movement between different localities [22, 23]. It has been shown that increased percentages of non-crop areas may lead to increased rates of parasitism and reduced rates of herbivory [24]. However, the spatial distribution of this non-crop vegetation, not only its percentage, is a relevant aspect for biological control, since parasitism rates can also be affected by the distance between local host populations [25]. Thus, appropriate landscape compositions may promote the balance between effective pest control and natural enemies conservation [26, 27].

In the attempt to understand how different landscape features influence species interactions and distributions, mathematical and computational models have long been recognised as important tools [28]. Considering more specifically host-parasitoid interactions, many studies were conducted to analyze stability [29–31], persistence [31] and spatial dynamics in parasitoid communities [32]. Although these studies have significantly contributed to our understanding of spatial aspects of host-parasitoid models, none of them investigated the effects of simultaneously considering different spatial scales on insect movement, nor its implications for the effectiveness of biological control in agroecosystems.

Here we investigate how the interplay between dispersal at different spatial scales and landscape structural connectivity and composition (in other words, physical configurations and spatial features) affect species dynamics and distributions. We propose a mathematical model to investigate the dynamics of host-parasitoid populations in a space characterized by two distinct scales of species dispersal. On the local scale, species interact on a landscape with regular structure and short distances between patches. On the regional scale, different landscapes are interconnected by dispersal through comparatively larger distances. Each of these landscapes is characterized by different distributions of non-crop habitats, where hosts are unable to grow and parasitoids are provided with alternative hosts. We call these habitats parasitoid refuges, directly influencing parasitoids’ spatial distributions and, consequently, their interaction dynamics with the primary hosts.

Specifically, we aim to assess (i) the direct impact of increasing proportions of parasitoid refuge areas on the host variables of density and occupancy, (ii) the differential effect of spatial distributions of these refuge areas on the host variables, and (iii) how regional connectivity patterns affect species dynamics. Additionally, we also evaluate the effect of the parasitoid dispersal radius on the host variables at different spatial scenarios.

Methods

Spatial structure

In our model, space is represented by two layers of networks, corresponding to two distinct scales, which we name external network and square lattice (or internal network) as we show in Fig 1.

[Figure omitted. See PDF.]

Fig 1. Local and regional spatial structures.

Illustration of spatial structuring in the model, showing two different scales. On the regional scale, connectivity is represented by an external network, where network nodes are connected as ring or star topologies. Each of these nodes corresponds to a regular square lattice on the local scale, with patches assigned two types: crop or refuge (non-crop) areas. For each proportion of refuge areas on the lattice, three refuge distributions are considered: random, block, and border.

https://doi.org/10.1371/journal.pone.0267037.g001

The external network represents the regional scale of an agricultural landscape, where each of the nodes corresponds to individual farms or production units, possibly linked. The links represent the possibility of dispersal of individuals between these units and, therefore, the network structure reflects the pattern of regional connectivity. We have used two different connectivity patterns for the external networks (Fig 1): (i) each node is linked only to its immediately adjacent neighbors (ring network), and (ii) several peripheral nodes are exclusively linked to one central node (star network). For the ring network, all the nodes are identical in relation to regional connectivity, whereas for the star network there is a clear separation in two groups (peripheral and central nodes).

The choice of these two arrangements highlights different spatial scenarios that might be imposed by regional topography: homogeneous connectivities, where all production units have the same number of neighbors, and heterogeneous connectivities, where the central unit behaves as a hub and the other units have only one link. The space between these units can represent urban areas, freshwater or any form of land use not suitable for the species we consider. While admittedly a simplification of reality, this representation enables the exploration of how differences in the structural connectivity of production units influence the dynamics of pests and natural enemies locally, emphasising the importance of integrated regional approaches to deal with common threats.

The internal networks describe the local structures of each of the production units and are presented as regular square lattices of size L x L. Each of the nodes on the square lattice represents a patch, formed by a small group of individual plants. The dispersal of individuals between different patches depends on the dispersal scale of each species (see description below).

Patches on the lattice can be of two types: crop or non-crop areas. We assume hosts to be specialists on that specific crop and, therefore, unable to grow on non-crop areas. Additionally, non-crop areas are also assumed to provide alternative hosts and food supplies for parasitoids, allowing reproduction and maintenance of small populations, and we refer to them as refuges. Note that this is different from the more common concept of host/prey refuge, where an organism obtains protection from predation by hiding. Here, the term refuge is used as a generalisation, referring to spatial features that allow increased survival of individuals, which, in our model, correspond only to parasitoids finding alternative hosts and resources at these locations. Henceforth, we use refuge as reference to parasitoid refuge in this context.

At the beginning of each simulation, a fixed fraction of the patches is assigned as refuge areas, while the rest remains as crop areas. We have chosen three spatial arrangements to model common patterns of non-crop vegetation in agricultural landscapes (Fig 1): (i) random, where a given fraction of refuge areas are randomly and homogeneously distributed on the lattice, (ii) block, where the refuge areas are grouped in a square cluster located at the centre of the lattice, and (iii) border, where refuge areas are randomly chosen at the perimeter of the lattice. The different types of patches, crop or refuge, present different local dynamics for each species, as described below.

Dynamics: Interaction

For the crop patches, interaction is given by a modified form of the host-parasitoid model proposed by Nicholson and Bailey [33], with parasitoids that attack hosts’ eggs and host density-dependent effects that become relevant at the larval stage. This formulation is relevant, for example, for the interaction between lepidopteran pests and their parasitoids, but not limited to them.

Male parasitoids are not explicitly included in the model, therefore the female qualifier will be omitted when referring to parasitoid populations. If Hi,t−1 and Fi,t−1 are the densities of hosts and parasitoids at patch i of the lattice at the previous generation (after dispersal step), the post-interaction host and parasitoid densities, hi,t and fi,t, will be given by:(1)where λ represents host growth rate, k its carrying capacity, and c the average number of parasitoid offspring emerging from each host. The term represents the fraction of hosts that escapes parasitoid attack and α gives the parasitoid attack rate.

Host dynamics can be written aswhere the function corresponds to the limitation of growth by density-dependent effects, analogous to the one initially proposed by Beverton and Holt [34]. The inclusion of the term on the function g(Hi,t−1, Fi,t−1), the main difference with the Berveton-Holt model, assumes an idea proposed by May [35]: that the number of hosts in the function for density dependence should be the one surviving parasitoid attack to the egg stage, therefore . The proposed formulation takes into account specific details of the interaction, without explicitly considering the life stages of the insects.

For the patches representing refuge areas, the interaction dynamics is given by:(2)

The subscript RA stands for refuge area. We assume hosts as specialist pests for the resources on crop areas. Therefore, host individuals that arrive at refuge site i, after a dispersal event on the previous generation, die due to the lack of resources. Parasitoids on the other hand are generalists and therefore able to survive and reproduce (with small growth rate λRA and carrying capacity kRA) by using alternative hosts.

If the densities fall below 1.0, they are set to 0.0 characterizing an event of local extinction. That site can be recolonized at a posterior step time.

Dynamics: Dispersal

Dispersal happens after the update to local populations due to the interaction. For each patch i a given fraction of hosts will disperse to neighboring patches. This fraction depends on the local density of conspecifics. Following Reigada et al. [36], the density of dispersing hosts leaving patch i at generation t is:(3)where μH is the maximum host dispersal rate in patches with large densities and h0 relates to the tolerance to conspecifics. For small values of h0, tends to μH hi,t, with dispersing hosts tending to the maximum fraction. For large values of h0, tends to 0 with no dispersing hosts per site.

Neighboring patches on the lattice can be reached by hosts with a dispersal radius RH centred on patch i. Besides, hosts can also disperse to adjacent lattices if the two lattices are connected on the external network. If dispersal happens to another lattice, a destination patch is randomly chosen on the border of the neighbor lattice for the hosts to arrive.

Dispersal is stochastic, with a subtraction of a unit of density to the density of dispersing hosts at each step. Each unit is sent to one of the patches inside the dispersal radius with probability , where is a normalization constant. Dispersal to neighboring lattices is given with probability , with RL representing the characteristic distance between nodes on the external network. The number of dispersal steps for each patch i is given by the integer part of (Eq 3). If no dispersal is implemented. After all patches i are accounted for on all lattices, host populations are updated.

Parasitoid dispersal is implemented in an analogous fashion, with a slight difference depending on the type of the local patch. For crop patches i, the density of dispersing parasitoids at generation t is given (according to Reigada et al. [36]) by:(4)but only if Hi,t > 0 (host density on site i after dispersal stage). If Hi,t = 0 all the local parasitoid population disperses, thus . Parameter μF gives the maximum parasitoid dispersal for small host densities and large parasitoid densities, while parameters H0 and f0 represent the necessary number of hosts to keep parasitoid at the patch, and the tolerance to parasitoid conspecifics, respectively.

For the refuge patch i, the density of dispersing parasitoids at generation t is:(5)

Parasitoids disperse within a radius RP, according to a similar process as the one for hosts, with neighboring patches being chosen with probability , where is a normalization constant. We assume stronger distance limitation on parasitoid dispersal, thus the squared distance on the denominator. Neighboring lattices can be reached with probability . Again, no dispersal step is implemented if . All parasitoid populations are updated after all sites i on all lattices are visited. Spatial population dynamics simulations were implemented in FORTRAN.

Parameterisation

We investigated two regional connectivity arrangements, ring and star networks, and three different distributions of refuge areas, random, block and border. A number N = 10 of nodes on the external network (regional connectivity) was considered for all the scenarios. For the square lattice, we fixed L = 33, resulting in a total L x L = 1089 patches for each lattice (if each patch has 10–100 m2, each square lattice has approximately 1–10 ha). We also assumed reflective boundary conditions for the square lattices. Six values of occupation fractions for refuge areas were considered in the lattices: 0.008, 0.023, 0.045, 0.074, 0.111, 0.155, and 0.207 (corresponding to 9, 25, 49, 81, 121, 169, and 225 patches), distributed according to one of the three arrangements. For the border distribution refuge area fractions were considered up to 0.111 (corresponding to an almost completely covered perimeter of the square lattice). Although all the analyses were done with the fractions of refuge areas, in the Results section (section) we chose to show them as percentages to improve clarity.

Parameters for the interaction and dispersal stages, in particular maximum dispersal rates and tolerances of hosts and parasitoids to conspecifics at the same patch (see complete description on Methods), were obtained from a previous host-parasitoid experimental setting [37]. Since we were interested in the long-term behaviour of the system, attack rate and initial densities were adjusted to guarantee species coexistence in the majority of the simulations. All simulations start with initial densities of 300 hosts and 4 parasitoids, placed in one of the crop patches. For the star network, individuals are initially placed at the central node of the star. For the ring network, the node index for the initial placement is irrelevant, since all nodes are identical in connectivity structure. Parasitoids have their dispersal radius parameterised according to short (RP = 1) or long (RP = 3) range dispersal. The set of parameters is summarized in Table 1.

[Figure omitted. See PDF.]

Table 1. Parameter values for the spatially-explicit simulation model.

https://doi.org/10.1371/journal.pone.0267037.t001

Variables of interest

The effects of landscape structure and parasitoid dispersal ranges on host dynamics were evaluated on two variables: the average host density per crop site (total host density on each lattice divided by the number of crop sites) and the average host occupancy per crop site (number of sites occupied by hosts divided by the number of crop sites on the lattice). Simulations were iterated for T = 5000 generations. Average host density and occupancy were calculated over the last 100 generations for each simulation. The value T = 5000 corresponds to a very conservative choice, since fewer than 100 initial steps are typically needed to reach average values of host densities at each node of the external network. This choice, however, ensures no initial transients affect the patterns we describe (see S1 Fig in S1 File for examples of typical time series of host and parasitoid densities). For each fraction of refuge areas considered, 50 sample simulations were run with the same set of parameters.

We also assessed how spatial structure and parasitoid traits influence time variation on host dynamics. We calculated the Coefficient of Variation (; ratio between the standard deviation, s, and the mean, , of a given sample) for the last 100 generations of the host density time series. We obtained one estimate of the CV for each of the 50 sample simulations, for each fraction of refuge areas and each value of parasitoid dispersal radius. The CV expresses how much the densities in a given interval deviate from the mean density for that interval, and can be used to compare scenarios with different means. Thus, we use it as a proxy to measure strong deviations from the mean density, or the propensity of a given scenario for the occurrence of population outbreaks.

Statistical analyses

In order to evaluate how the two dependent variables of interest, namely average host density per crop site and average host occupancy per crop site, would respond to changes in percentage and distribution of refuge areas, as well as to the network structure representing regional connectivity, we have performed linear regression analyses to evaluate significant differences in the slopes. The fraction of refuge areas was used as continuous independent variable with the three distributions of refuge areas (random, block, and border) as levels of a categorical variable. For the star network, an additional categorical variable representing the connectivity degree of nodes (peripheral and central) was also considered. Comparisons between slopes obtained from linear regressions were also performed for the coefficient of variation of host density.

For the comparisons, we have performed T-tests for statistical significance of the difference in the slopes obtained from the regressions, considering each parasitoid radius (RP = 1 and RP = 3) separately. Thus, for the ring network we have pairwise comparisons for the slopes for each categorical distribution of refuge areas, while for the star network, pairwise comparisons for the distributions are done for each category of connectivity separately. Homoscedasticity and normality of the residuals for the regressions were checked for consistency. All the statistical analyses were performed in R and the results can be found in the Supporting Information (see S1 and S2 Tables in S1 File).

Results

Ring network

Fig 2 shows the results for the average host density per crop site (Fig 2—first row) and the average host occupancy per crop site (Fig 2—second row) for the ring network, as functions of the percentage of refuge areas on the square lattices. For each percentage of refuge areas, results were obtained for two parasitoid dispersal radii RP = 1 and RP = 3. Ensemble averages are shown with error bars, corresponding to one standard error of the mean. For each set of parameters, 50 samples were simulated. Densities and occupancies for the 10 nodes of the external network are shown, with 10 circles (RP = 1) and 10 triangles (RP = 3) for each value of percentage of refuge areas.

[Figure omitted. See PDF.]

Fig 2. Host variables—Ring network.

Average host density per crop site (first row) and average host occupancy per crop site (second row) as functions of the percentage of refuge areas on the lattice, for random, block, and border distributions of refuge areas (first, second, and third columns, respectively). Ten points are shown for each percentage of refuge areas, one for each node on the external network. Each point corresponds to the average density or occupancy for that node over 50 simulations, with error bars showing one standard error of the mean. Results are shown for two values of parasitoid dispersal radius RP = 1 (blue circles) and RP = 3 (orange triangles). External network (regional connectivity) has ring topology.

https://doi.org/10.1371/journal.pone.0267037.g002

For all the scenarios, the 10 points for each fraction of refuge areas present very similar values, showing no significant difference in host variables among different nodes of the external network. This result follows directly from the topological equivalence between nodes on the ring network, with enough time given to populations to fully spread through the network (averages were calculated for the last 100 generations of the time series). Values for density and occupancy of hosts for RP = 3 are always lower than the values obtained for RP = 1.

For all the distributions of refuge areas, the average density of hosts per crop site decreases as the fraction of refuges increases, with the largest reduction obtained for the random distribution (Fig 2, top left panel). The decrease in host density for the random distribution, however, is followed by the largest increase in average occupancy of hosts per crop site on the lattice (Fig 2, bottom left panel). This increase does not represent a simple spread of host individuals on the lattice, since the total lattice densities (host density per site times number of crop sites) decrease when the fraction of refuge areas goes from 0.008 to 0.207.

Each distribution of refuge areas determines different responses of the dependent variables as the proportion of refuge areas increases. We show significant differences on the slopes of the regression curves for both host density and host occupancy, comparing the three refuge distributions in pairs (S1 Table in S1 File), for both values of parasitoid dispersal ranges (RP = 1 and RP = 3). The only exception are the slopes for host density of block and border distributions for RP = 1 (Fig 2, top central and top right panels) which are not significantly different (p = 0.4043—S1 Table in S1 File).

Star network

Fig 3 shows the effect of the distributions of refuge areas on host density and occupancy as refuge area increases, for both parasitoid dispersal ranges (RP = 1 and RP = 3) for the star network. Host density and host occupancy show a clear separation of the values in two distinct categories: the peripheral nodes of the external network (empty symbols) and the central node (filled symbols). The central node presents consistently lower densities and higher occupancies compared to peripheral nodes, for each percentage of refuge areas on all the scenarios.

[Figure omitted. See PDF.]

Fig 3. Host variables—Star network.

Average host density per crop site (first row) and average host occupancy per crop site (second row) as functions of the percentage of refuge areas on the lattice, for random, block, and border distributions of refuge areas (first, second, and third columns, respectively). Ten points are shown for each percentage of refuge areas, one for each node on the external network. Peripheral nodes are represented as empty symbols and central nodes as filled symbols. Each point corresponds to the average density or occupancy for that node over 50 simulations, with error bars showing one standard error of the mean. Results are shown for two values of parasitoid dispersal radius RP = 1 (blue circles) and RP = 3 (orange triangles). External network (regional connectivity) has star topology.

https://doi.org/10.1371/journal.pone.0267037.g003

A reduction of the difference between peripheral and central nodes as the fraction of refuge areas increases can be seen, in particular for the random and border distributions (both for density—Fig 3, top left and top right panels—and occupancy—Fig 3, bottom left and bottom right panels). For the block distribution, this reduction in the differences between node categories is less pronounced and shows a slight increase for host occupancy at RP = 3 (Fig 3, bottom central panel).

The effects of refuge distributions on the responses of host densities and host occupancies to the percentage of refuge areas are statistically significant, separately compared for each parasitoid radius and node category (S2 Table in S1 File). The exceptions are the difference between block and border distributions on host density, for central nodes and RP = 1, and the difference between random and border on host occupancy, for peripheral nodes and RP = 3, which present large confidence intervals on the estimation of slopes.

The separation between central and peripheral nodes reveals an important aspect that could not be observed on the ring network. For the ring network, the border distribution presents a small variation in host density and occupancy, as the fraction of refuge areas increases (Fig 2, top right and bottom right panels). For the same distribution of refuge areas on the star network, there is a clear decrease in density and increase in occupancy for the peripheral nodes as the percentage of refuge areas increases, while the central nodes present no clear variation (Fig 3, top right and bottom right panels). Thus, the connectivity of this lattice on the regional scale influences how host densities and occupancies will be affected by local percentages of refuge areas. For other refuge distributions (random and block) in star networks, host densities decrease and host occupancies increase as refuge areas increase, for peripheral and central nodes, with significant differences in these tendencies when refuge distributions are compared (S2 Table in S1 File).

Values for host density and occupancy corresponding to the parasitoid dispersal radius RP = 3 are always lower than the ones obtained for RP = 1 (Fig 3), if comparisons are made between nodes on the same category of regional connectivity (peripheral nodes with RP = 1 compared to peripheral nodes with RP = 3, for example).

Coefficients of variation (CV) of host density

Fig 4 shows the results for the ensemble averages of the CV for host density per crop site, for random, block and border distributions of refuge areas, in ring and star external networks. Our results show that the interplay between regional connectivity and spatial distribution of refuges influences the variability in relation to the mean density (measured by the CV) for the majority of the scenarios, with significant differences in regression slopes when comparing pairs of distributions of refuge areas (S1 and S2 Tables in S1 File). Important exceptions are the comparisons between block and border, for the ring network with RP = 3, and between block and border in the star network, for both central and peripheral nodes, with RP = 1, for which the differences in slopes are not significant.

[Figure omitted. See PDF.]

Fig 4. Coefficients of variation—Ring and star networks.

Coefficients of variation (CV) for host density per crop site as a function of the percentage of refuge areas for the ring (first row) and star (second row) networks with random, block, and border distributions of refuge areas (first, second, and third columns, respectively). Ten points are shown for each percentage of refuge areas, one for each node on the external network. For the star network, peripheral nodes are represented as empty symbols and central nodes as filled symbols. Each point corresponds to the average CV for that node over 50 simulations, with error bars representing one standard error of the mean. Results are shown for two values of parasitoid dispersal radius RP = 1 (blue circles) e RP = 3 (orange triangles).

https://doi.org/10.1371/journal.pone.0267037.g004

For the random distribution, an increase in the percentage of refuges leads to a decrease in the CV for both ring and star networks. The variation on the CV is not as evident for block and border distributions. For the star network, however, a clear separation in the values of CV for peripheral and central nodes (empty and filled symbols, respectively) is also shown, with central nodes presenting consistently lower values of the CV, for every percentage of refuge areas. For the border distribution in the star network (Fig 4, bottom right panel), an initial decrease in CV can also be noted at peripheral nodes, stabilising after 5% of refuge areas. For the central nodes, the CV increases following an increase in refuge areas, leading to a reduction on the separation between peripheral and central nodes.

The values of the CV for RP = 3 are slightly greater than the ones for RP = 1 in all of the scenarios, except for the border distribution in the ring network (Fig 4, top right panel), and for larger fractions of refuge areas for the border distribution in the star network (Fig 4, bottom right panel).

Discussion

Understanding how the structural characteristics of the landscapes affect species dynamics and distribution is fundamental to devise better strategies for pest control and conservation. Landscape composition has been widely recognised to have an important effect on biological control in agroecosystems, specifically in terms of proportion and distribution of non-crop or semi-natural habitats, by providing shelter and alternative resources for natural or introduced enemies [14, 22, 38, 39]. Moreover, it has also been shown that different scales of landscape structure from local to regional level are important to understand pest dynamics and distributions in natural and managed environments [7, 8, 40]. Our main goal was to establish how these two scales of spatial organisation interact and determine different responses to population dynamics of interacting species.

We used a host-parasitoid mathematical model to investigate the influence of spatial structure on population dynamics, focusing primarily on host density and occupancy on the landscape. We explored two different spatial scales: at the local scale, analysing different proportions and distribution of refuge areas, and at the regional scale, different patterns of connectivity on the landscape. Parasitoid dispersal radius, another relevant parameter for the host-parasitoid interaction, was also varied, showing how the interplay between biotic and abiotic aspects might produce both qualitative and quantitative changes on the observed patterns. Our results reveal three non-intuitive key findings about the role played by refuge areas on species dynamics.

First, we show that hosts’ temporal and spatial dynamics are dependent not only on the proportion of parasitoid refuges on the landscape, but also on their distribution and on the connectivity structure at the regional scale. In our model, refuge areas correspond to non-crop patches where it is possible to find alternative hosts and food supplies, allowing reproduction and maintenance of a small population of parasitoids [22]. Therefore, it was expected that an increase in the number of refuge sites (viewed as landscape sources of parasitoids) would always be associated with decreases in density and occupancy of hosts at crop sites [41]. However, our results show that decreases in density and occupancy are not always the rule, with significant differences on the response of these variables to refuge proportions, depending on the spatial distributions of refuges.

Second, the results we have obtained show how host populational patterns are dependent on the parasitoid dispersal radius, but also how spatial structure plays a fundamental role. While in most cases an increased parasitoid dispersal radius leads to decreased values for host density and occupancy and increased values for coefficients of variation, changes to this general expectation can be seen for specific spatial arrangements (see border distribution in Fig 4, top right and bottom right panels). This fact suggests that the choice of agents for biological control should take into account not only dispersal or interaction traits, but also the landscape structure where this control is meant to be applied [42].

Third, the difference on the patterns observed for ring and star networks shows that the same distribution of refuge areas at the local scale may present meaningful differences on the dynamical effects depending on the landscape connectivity at the regional scale. The separation on the host population variables related to peripheral and central nodes for the star network highlights an important interplay between local distribution of refuge areas and landscape connectivity, since the more connected nodes in the regional scale (central nodes on the external networks) consistently present smaller values for host density and larger values for host occupancy. For the ring network, there was no significant distinction between patterns observed for each node.

In an agricultural context, where hosts could represent insect-pests for crop areas, our results might inform design strategies for biological control, through the analysis of landscape connectivity and restructuring of refuges to sustain natural or introduced populations. With units of production located on a regional landscape with different degrees of connectivity, biological control of pests would benefit of an integrated regional approach [43, 44]. In this scenario, our results suggest that the identification of central and peripheral nodes, might help to articulate directed actions on key units, such as engineering refuge distributions or augmented biological control.

Comparing the three refuge distributions, the random distribution presented the largest reduction on host density as a function of the proportion of refuge areas, in agreement to the beneficial effect of high edge density on pest control in real landscapes [10]. The choice of this refuge distribution as the best candidate for pest control in the context of our modelling approach, however, must be treated with caution. As a simplifying assumption, we have not considered possible variations of population growth of parasitoids at refuge sites. At these sites, we assume logistic growth for parasitoids, with population stochasticity only due to dispersal. A complex population model could be obtained, for example, by considering the dynamics of alternative hosts at refuges. This would have a relevant effect on the dynamical response of populations subject to the random distribution of non-crop areas, particularly for small fractions of refuges on the lattice, since these sites would be fairly isolated. In this scenario, the assumption that all of them could constitute, at any time, refuges capable of supporting stable communities of alternative hosts would no longer hold. Thus, the dynamics of alternative hosts at refuge areas constitutes an important aspect to be addressed in future work.

We show a clear distinction in effectiveness of host density reduction, measured by differences in slopes for different spatial distributions of parasitoid refuges. This relates to a long running discussion regarding the relevance of refuge and native area fragments on the maintenance of natural enemies populations and biodiversity conservation in general [45–48], specially regarding the effects of size and distribution of these fragments. Although our work hints to the relative role of fragmentation for biological control on agroecosystems, a complete evaluation of the matter would also consider the dispersal patterns of the species under study [49]. Nevertheless, our results clearly present a separation between the percentage of refuge habitats and their spatial distribution, which has also been the focus of much recent debate [50, 51].

Finally, although we have focused essentially on the population variables of hosts, understanding the reciprocal actions of these variables with the population dynamics of parasitoids is extremely important to comprehend the biotic and abiotic factors that condition the spatio-temporal species patterns in this context. In particular, future work should address how refuge distributions and regional connectivity influence the occurrence of population outbreaks for both species in time and space. The study of these interrelations can provide valuable knowledge to inform policy on pest control and spatial structuring of these landscapes.

Supporting information

S1 File. Supporting information with time series for host and parasitoid densities, for different distributions of refuge areas (S1 Fig) and summary statistics for ring and star networks (S1 & S2 Tables).

https://doi.org/10.1371/journal.pone.0267037.s001

(PDF)

Citation: D. Fernandes L, Mata AS, Godoy WAC, Reigada C (2022) Refuge distributions and landscape connectivity affect host-parasitoid dynamics: Motivations for biological control in agroecosystems. PLoS ONE 17(4): e0267037. https://doi.org/10.1371/journal.pone.0267037

References

1. Lustig A, Stouffer DB, Doscher C, Worner SP. Landscape metrics as a framework to measure the effect of landscape structure on the spread of invasive insect species. Landscape Ecology. 2017;32(12):2311–2325.

2. Thomas C. Extinction, colonization, and metapopulations: environmental tracking by rare species. Conservation Biology. 1994;8(2):373–378.

3. Sõber V, Leps M, Kaasik A, Mänd M, Teder T. Forest proximity supports bumblebee species richness and abundance in hemi-boreal agricultural landscape. Agriculture, Ecosystems & Environment. 2020;298:106961.

4. Urban D, Keitt T. Landscape connectivity: a graph-theoretic perspective. Ecology. 2001;82(5):1205–1218.

5. Fortuna MA, Gómez-Rodríguez C, Bascompte J. Spatial network structure and amphibian persistence in stochastic environments. Proceedings of the Royal Society B: Biological Sciences. 2006;273(1592):1429–1434. pmid:16777733

6. Altermatt F, Bieger A, Carrara F, Rinaldo A, Holyoak M. Effects of connectivity and recurrent local disturbances on community structure and population density in experimental metacommunities. PLoS One. 2011;6(4):e19525. pmid:21559336

7. Seidl R, Müller J, Hothorn T, Bässler C, Heurich M, Kautz M. Small beetle, large-scale drivers: how regional and landscape factors affect outbreaks of the European spruce bark beetle. Journal of Applied Ecology. 2016;53(2):530–540.

8. Senf C, Campbell EM, Pflugmacher D, Wulder MA, Hostert P. A multi-scale analysis of western spruce budworm outbreak dynamics. Landscape ecology. 2017;32(3):501–514.

9. Raffa KF, Aukema BH, Bentz BJ, Carroll AL, Hicke JA, Turner MG, et al. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: the dynamics of bark beetle eruptions. Bioscience. 2008;58(6):501–517.

10. Martin EA, Dainese M, Clough Y, Báldi A, Bommarco R, Gagic V, et al. The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe. Ecology letters. 2019;22(7):1083–1094. pmid:30957401

11. Karp DS, Chaplin-Kramer R, Meehan TD, Martin EA, DeClerck F, Grab H, et al. Crop pests and predators exhibit inconsistent responses to surrounding landscape composition. Proceedings of the National Academy of Sciences. 2018;115(33):E7863–E7870. pmid:30072434

12. Medeiros HR, Martello F, Almeida EA, Mengual X, Harper KA, Grandinete YC, et al. Landscape structure shapes the diversity of beneficial insects in coffee producing landscapes. Biological Conservation. 2019;238:108193.

13. Malard JJ, Adamowski JF, Díaz MR, Nassar JB, Anandaraja N, Tuy H, et al. Agroecological food web modelling to evaluate and design organic and conventional agricultural systems. Ecological Modelling. 2020;421:108961.

14. Le Gal A, Robert C, Accatino F, Claessen D, Lecomte J. Modelling the interactions between landscape structure and spatio-temporal dynamics of pest natural enemies: Implications for conservation biological control. Ecological modelling. 2020;420(15):108912.

15. Wajnberg E, Bernstein C, Van Alphen J. Behavioral ecology of insect parasitoids: from theoretical approaches to field applications. John Wiley & Sons; 2008.

16. Bianchi FJ, Booij C, Tscharntke T. Sustainable pest regulation in agricultural landscapes: a review on landscape composition, biodiversity and natural pest control. Proceedings of the Royal Society B: Biological Sciences. 2006;273(1595):1715–1727. pmid:16790403

17. Veres A, Petit S, Conord C, Lavigne C. Does landscape composition affect pest abundance and their control by natural enemies? A review. Agriculture, Ecosystems & Environment. 2013;166:110–117.

18. González E, Štrobl M, Janšta P, Hovorka T, Kadlec T, Knapp M. Artificial temporary non-crop habitats support parasitoids on arable land. Biological Conservation. 2022;265:109409.

19. Barbosa P, Benrey B. The influence of plants on insect parasitoids: implications for conservation biological control. In: Conservation biological control. Elsevier; 1998. p. 55–82.

20. Bellone D, Björkman C, Klapiwijk M. Top-down pressure by generalist and specialist natural enemies in relation to habitat heterogeneity and resource availability. Basic and Applied Ecology. 2019;43.

21. Letourneau DK, Armbrecht I, Rivera BS, Lerma JM, Carmona EJ, Daza MC, et al. Does plant diversity benefit agroecosystems? A synthetic review. Ecological applications. 2011;21(1):9–21. pmid:21516884

22. Rusch A, Valantin-Morison M, Sarthou JP, Roger-Estrade J. Biological control of insect pests in agroecosystems: effects of crop management, farming systems, and seminatural habitats at the landscape scale: a review. In: Advances in agronomy. vol. 109. Elsevier; 2010. p. 219–259. https://doi.org/10.1016/B978-0-12-385040-9.00006-2

23. Fraser SE, Dytham C, Mayhew PJ. Determinants of parasitoid abundance and diversity in woodland habitats. Journal of applied ecology. 2007;44(2):352–361.

24. Schmidt MH, Thies C, Tscharntke T. Landscape context of arthropod biological control. In: Gurr GM, Wratten SD, Altieri MA, editors. Ecological engineering for pest management: advances in habitat manipulation for arthropods. CSIRO Publishing; 2004. p. 55–63.

25. Pareja M, Brown VK, Powell W. Aggregation of parasitism risk in an aphid-parasitoid system: Effects of plant patch size and aphid density. Basic and Applied Ecology. 2008;9(6):701–708.

26. Roderick GK, Navajas M. Genes in new environments: genetics and evolution in biological control. Nature Reviews Genetics. 2003;4(11):889. pmid:14634636

27. Souza B, Vazquez Luiz L, Marucci Rosangela C. Natural Enemies of Insect Pests in Neotropical Agroecosystems. Springer; 2019.

28. Newman EA, Kennedy MC, Falk DA, McKenzie D. Scaling and complexity in landscape ecology. Frontiers in Ecology and Evolution. 2019;7:293.

29. Murdoch WW, Briggs CJ, Nisbet RM, Gurney WS, Stewart-Oaten A. Aggregation and stability in metapopulation models. The American Naturalist. 1992;140(1):41–58. pmid:19426065

30. Briggs CJ, Hoopes MF. Stabilizing effects in spatial parasitoid–host and predator–prey models: a review. Theoretical population biology. 2004;65(3):299–315. pmid:15139366

31. Vandermeer J, Jackson D. Stabilizing intransitive loops: self-organized spatial structure and disjoint time frames in the coffee agroecosystem. Ecosphere. 2018;9(12):e02489.

32. Amarasekare P. Spatial dynamics in a host–multiparasitoid community. Journal of Animal Ecology. 2000;69(2):201–213.

33. Nicholson AJ, Bailey VA. The Balance of Animal Populations.—Part I. Journal of Zoology. 1935;105(3):551–598.

34. Beverton RJ, Holt SJ. On the dynamics of exploited fish populations, Fishery Investigations Series II, Vol. XIX, Ministry of Agriculture. Fisheries and Food. 1957;1:957.

35. May R, Hassell M, Anderson R, Tonkyn D. Density dependence in host-parasitoid models. The Journal of Animal Ecology. 1981; p. 855–865.

36. Reigada C, de Aguiar MAM, Fernandes LD. Demographic Processes in Spatially Structured Host-Parasitoid Systems. In: Ecological Modelling Applied to Entomology. Springer; 2014. p. 11–38.

37. Reigada C, de Aguiar MA. Host–parasitoid persistence over variable spatio-temporally susceptible habitats: bottom–up effects of ephemeral resources. Oikos. 2012;121(10):1665–1679.

38. Perez-Alvarez R, Nault BA, Poveda K. Effectiveness of augmentative biological control depends on landscape context. Scientific reports. 2019;9(1):1–15. pmid:31209256

39. Chaplin-Kramer R, Kremen C. Pest control experiments show benefits of complexity at landscape and local scales. Ecological Applications. 2012;22(7):1936–1948. pmid:23210310

40. Wildemeersch M, Franklin O, Seidl R, Rogelj J, Moorthy I, Thurner S. Modelling the multi-scaled nature of pest outbreaks. Ecological Modelling. 2019;409:108745.

41. Zhao ZH, He DH, Hui C. From the inverse density–area relationship to the minimum patch size of a host–parasitoid system. Ecological research. 2012;27(2):303–309.

42. Haan NL, Zhang Y, Landis DA. Predicting landscape configuration effects on agricultural pest suppression. Trends in Ecology & Evolution. 2020;35(2):175–186. pmid:31699410

43. Tscharntke T, Bommarco R, Clough Y, Crist T, Kleijn D, Rand T, et al. Conservation biological control and enemy diversity on a landscape scale [Erratum: 2008 May, v. 45, issue 2, p. 238-253.]. Biological control: theory and application in pest management. 2007;.

44. Thies C, Tscharntke T. Landscape structure and biological control in agroecosystems. Science. 1999;285(5429):893–895. pmid:10436158

45. Quinn JF, Harrison SP. Effects of habitat fragmentation and isolation on species richness: evidence from biogeographic patterns. Oecologia. 1988;75(1):132–140. pmid:28311846

46. Andren H. Effects of habitat fragmentation on birds and mammals in landscapes with different proportions of suitable habitat: a review. Oikos. 1994; p. 355–366.

47. Tscharntke T, Steffan-Dewenter I, Kruess A, Thies C. Characteristics of insect populations on habitat fragments: a mini review. Ecological research. 2002;17(2):229–239.

48. Fahrig L. Ecological responses to habitat fragmentation per se. Annual Review of Ecology, Evolution, and Systematics. 2017;48:1–23.

49. Thompson SE, Chisholm RA, Rosindell J. Characterising extinction debt following habitat fragmentation using neutral theory. Ecology Letters. 2019;22(12):2087–2096. pmid:31612627

50. Fletcher RJ Jr, Didham RK, Banks-Leite C, Barlow J, Ewers RM, Rosindell J, et al. Is habitat fragmentation good for biodiversity? Biological conservation. 2018;226:9–15.

51. Fahrig L, Arroyo-Rodríguez V, Bennett JR, Boucher-Lalonde V, Cazetta E, Currie DJ, et al. Is habitat fragmentation bad for biodiversity? Biological Conservation. 2019;230:179–186.