Combining DNA metabarcoding and ecological networks to inform conservation biocontrol by small vertebrate predators.

In multifunctional landscapes, diverse communities of flying vertebrate predators provide vital services of insect pest control. In such landscapes, conservation biocontrol should benefit service-providing species to enhance the flow, stability and resilience of pest control services supporting the production of food and fiber. However, this would require identifying key service providers, which may be challenging when multiple predators interact with multiple pests. Here we provide a framework to identify the functional role of individual species to pest control in multifunctional landscapes. First, we used DNA metabarcoding to provide detailed data on pest species predation by diverse predator communities. Then, these data were fed into an extensive network analysis, in which information relevant for conservation biocontrol is gained from parameters describing network structure (e.g., modularity) and species roles in such network (e.g., centrality, specialization). We applied our framework to a Mediterranean landscape, where 19 bat species were found to feed on 132 insect pest species. Metabarcoding data revealed potentially important bats that consumed insect pest species in high frequency and/or diversity. Network analysis showed a modular structure, indicating sets of bat species that are required to regulate specific sets of insect pests. A few generalist bats had particularly important roles, either at network or module levels. Extinction simulations highlighted six bats, including species of conservation concern, which were sufficient to ensure that over three-quarters of the pest species had at least one bat predator. Combining DNA metabarcoding and ecological network analysis provides a valuable framework to identify individual species within diverse predator communities that might have a disproportionate contribution to pest control services in multifunctional landscapes. These species can be regarded as candidate targets for conservation biocontrol, although additional information is needed to evaluate their actual effectiveness in pest regulation.

Introduction

Landscapes providing multiple products and services (i.e., multifunctional) are key to combine food security and biodiversity conservation worldwide, as they integrate human activities with the preservation of ecosystem structure and function (Kremen and Merenlender 2018, Manning et al. 2018). In such landscapes, biodiversity contributes to food and fiber production, and to reduce its negative environmental externalities, by providing critical services such as pollination and pest control (Mace et al. 2012, Maas et al. 2016, Kaiser‐Bunbury et al. 2017, Heath and Long 2019). The management of multifunctional landscapes should therefore target at benefiting service‐providing organisms to enhance the flow, stability and resilience of ecosystem services supporting production (Doré et al. 2011, Bommarco et al. 2013, Tittonell 2014). An important example is conservation biocontrol, which aims at increasing populations of natural enemies to reduce crop and forestry losses by pests (Shields et al. 2019).

Designing conservation biocontrol strategies is particularly challenging in multifunctional landscapes, because multiple land uses coexist across space and over time (e.g., crop rotation), with each land use associated with multiple pests, and each pest predated by multiple predators (Bianchi et al. 2006, Chaplin‐Kramer et al. 2011). Further complexity is added by highly mobile predators such as small flying vertebrates, which can provide vital biocontrol services while feeding across land use types on multiple pest species (Maine and Boyles 2015, Maas et al. 2016). In these circumstances, regulation of all pests represented in a multifunctional landscape may require diverse predator communities rather than any particular predator species, and generic management prescriptions such as maintaining landscape heterogeneity and/or maximizing coverage by natural habitats may be sufficient to enhance pest biocontrol (Bianchi et al. 2006, Rusch et al. 2016, Karp et al. 2018). However, targeted efforts to benefit individual species may still be needed, because some particular species may have a disproportionate importance for pest control (Cleveland et al. 2006, Wanger et al. 2014). This is particularly relevant, for instance, for habitat and diet specialists, and of endangered or otherwise vulnerable species, which may be at risk of declining and therefore of releasing specific pests from effective biocontrol. Understanding the role of individual species for pest suppression in multifunctional landscapes is therefore essential to inform conservation biocontrol.

Species' roles in systems involving multiple predators, multiple pests and multiple land use types are challenging to evaluate. For instance, while exclosure experiments have been used successfully to show the importance of bats and birds to pest suppression in individual crops (Maine and Boyles 2015, Maas et al. 2016), the approach is less suited for identifying the actual species involved, if not combined with other techniques, and it may be hard to apply simultaneously to multiple crops in fragmented and diverse landscapes. A practical alternative to approximate this problem is the analysis of species interaction networks, which can aid the identification of key species needed for community functioning (Harvey et al. 2017), and therefore help to identify targets for conservation biocontrol (Table 1). In bipartite networks such as those representing predator‐pest interactions, a species having many interactions (i.e., high degree centrality) can be considered particularly important, especially if the network has a nested structure (i.e., the interactions of species with a lower degree are a subset of those with a higher degree; Bascompte et al. 2003). However, if the network is modular (i.e., with groups of species interacting strongly with other species within but not across modules), species acting as “network hubs” (i.e., interacting with many species across modules) may be the most important, although relevance should also be given to “module hubs” (i.e., interacting with many species within a module; Delmas et al. 2019). Attention by managers should also be given to species that have low trophic niche overlap with others and therefore complement their functional role, although functionally redundant species may also contribute to the stability of pest control services (Biggs et al. 2020). Therefore, although the analysis of ecological network cannot provide definite information on the actual pest suppression role of individual species, it may provide important clues on candidate species fulfilling such role.

Table 1

Conceptual framework for combining DNA metabarcoding and the analysis of predator‐pest networks to inform conservation biocontrol in multifunctional landscapes.

Metric	Concept	Implications
Diet
Species richness (R: 0‐n)	Number of different pest species consumed by a predator	Predators with high pest richness (generalists) may be effective at regulating a wide range of pests
Frequency of interaction (FI: 0‐1)	Proportion of interactions with a pest species in relation to the total number of interactions with all prey species	Predators with high FI include a large proportion of pests in the diet and may be particularly relevant for pest regulation
Frequency of occurrence (FO: 0‐1)	Proportion of individuals consuming a pest species in relation to the total number of individuals analyzed	Predators with high FO include a large number of individuals consuming pest species and may be particularly relevant for pest regulation
Diet distinctiveness (DD: 0‐1)	Level of distinctiveness of a predator' diet, computed as one minus the average of pairwise Pianka's niche overlaps between the focal and each other species	Predators with high diet distinctiveness may regulate pests that are covered by no other species
Network structure
Modularity1 (Q: 0‐1)	Modules identify aggregated sets of interacting predators and pests, when within‐module interactions are more prevalent than between‐module interactions	Predators from different modules need to be represented in a landscape to control pests from their corresponding modules
Nestedness2 (wNODF: 0‐1)	Nestedness indicates the extent to which specialist predators interact with proper subsets of the pest species interacting with generalists	In nested networks the pest control services provided by specialists are redundant to those provided by generalists
Network specialization3 (H 2′: 0‐1)	Indicates the extent to which a network is dominated by specialist predators (i.e., number of interactions lower than expected by chance)	In specialized networks, a large number of complementary predator species may be needed to provide biocontrol for all pest species
Species roles in the network
Centrality
Normalized Degree4 (ND: 0‐1)	Indicates the number of interactions per predator species (degree) divided by the number of possible interacting pest species	Predators with high ND are generalists that may be effective at regulating a wide range of pests, and that are core in predator‐pest network structure and enhance its robustness
Betweenness4 (BC: 0‐1)	BC measures the importance of a node as a connector between different parts of the network	Predators with high BC may be important to regulate pest species that otherwise are consumed by few other predators
Closeness4 (CC: 0‐1)	CC measures the proximity of a species to all other species in the network	Predators with high CC share many interaction partners and may be important to provide redundancy to predator‐pest interactions
Module‐based roles
Within‐module degree4,5 (z: ‐n‐n)	Measures within‐module connectivity, and therefore the importance of a species to its own module; species with z > a threshold computed from null models are considered “hubs” and have significantly more interactions within their module than expected by chance	Predators with high z may be important because they potentially regulate a large range of species within its module
Participation coefficient4,5 (c: 0‐1)	Measures between‐modules connectivity, and therefore the distribution of focal species' interactions across modules; species with c > a threshold computed from null models are considered “connectors” and have significantly more interactions across modules than expected by chance	Predators with high c may be important because they contribute to the cohesion of the network, being able to regulate pest species from different modules
Module‐based species classifications4,5	(1) Module hubs – high z and low c; (2) network hubs – high z and c, assuring cohesion of the network and cohesion of a module; (3) low z and c – peripherals; low z and high c – connectors that “glue” different modules together	Predators that are network hubs may be particularly important to regulate a wide range of pest species
Specialization
Standardized specialization index6 (d′: 0‐1)	Measures specialization as discrimination from expectation based on how many interactions a predator has. High d′ values indicate species of a trophic level that interact with few species and with species that few (or no) other species interact	Specialist predators with high d′ may be particularly important because they are less redundant than generalists, and therefore may have unreplaceable roles at regulating some pest species
Vulnerability to extinction
Simulated sequential extinction curves7	Indicates how the loss of a predator species can release pest species from biocontrol by the studied community, based on simulated extinctions following a predefined sequence (e.g., random, high‐ to low‐linked species, etc.) Takes the strong assumptions that prior to extinction the predator exerted a top‐down control on their pest prey species; that the predator needs to go extinct for the pest to be released; and that there is no “rewiring” of the network (i.e., predators taking each other's place in response to extinctions)	Predators that greatly affect the shape of the extinction curve can be critical to assure pest control in the community

For each metric we provide in brackets its acronym and range of variation.

1Beckett (2016); 2Almeida‐Neto and Ulrich (2011); Dormann et al. (2009); 3 Blüthgen et al. (2006); 4Cirtwill et al. (2018); 5Hackett et al. (2019); 6Dormann (2011); 7Memmott et al. (2004).

Network analysis has been widely used to understand species roles in mutualistic (e.g., pollination and seed dispersal) and antagonistic (e.g., parasitoid‐host and herbivory) networks, but much less attention has been given to predator‐pest interaction networks, due to difficulties in assessing the diet of natural predator communities (but see Roubinet et al. 2018, Feit et al. 2019, Sint et al. 2019). Molecular diet analysis through DNA metabarcoding has recently overcome this problem (Pompanon et al. 2012, Nielsen et al. 2018), by simultaneously providing species level identification of hundreds of prey and the capacity to processing several hundred samples in a short period, and therefore allowing the reconstruction of insectivore diets with relative ease (Galan et al. 2018, Gordon et al. 2019). Moreover, reduced costs, coupled with technical refinements such as multimarker approaches (da Silva et al. 2019) and the availability of ever more comprehensive barcode databases (e.g., Dincă et al. 2015, Hendrich et al. 2015), are making this method increasingly powerful and readily available for practical applications. Using DNA metabarcoding and ecological networks to inform landscape and conservation management has already been proposed (Evans et al. 2016, Clare et al. 2019), but the approach remains largely unexplored.

This study provides a framework combining DNA metabarcoding and ecological network analysis to identify the functional role of individual species to pest control in multifunctional landscapes. We illustrate its practical application with a case study focusing on a diverse bat community inhabiting a Mediterranean landscape comprising multiple crops and forest production systems, which can be attacked by over one hundred insect pests. We focused on bats because they play important roles as arthropod pest suppressors (Maine and Boyles 2015, Maas et al. 2016), while Mediterranean mosaic landscapes are often considered the epitome of multifunctionality (Bugalho et al. 2011). Specifically, the study aimed at: (1) describing predation on arthropod pests (i.e., frequency of occurrence and interaction) by all bat species; (2) using this information to produce a predator‐pest network and characterize its properties; (3) estimating species' roles in the network by assessing species' centrality and module inter‐ and intraconnectivity; and (4) estimating the vulnerability of the service to bat extinctions, in particular of species of conservation concern. Results highlight the potential of DNA metabarcoding combined with ecological network analysis to inform conservation biocontrol in systems involving multiple predators and pests.

Materials and Methods

Study area

Bats were sampled in north‐eastern Portugal, within the Baixo Sabor Long Term Ecological Research Site (https://deims.org/45722713‐80e3‐4387‐a47b‐82c97a6ef62b), encompassing the watershed of the lower reaches of the Sabor River and the catchment of its tributaries. Elevation ranges between 100 m and 1,100 m asl. The area is included in the Meso‐Mediterranean and Supra‐Mediterranean zones, with 500–1,000 mm of annual rainfall, and 10–16°C of mean annual temperature. Land cover is dominated (>80%) by Mediterranean oak forests (Quercus spp.) and shrublands (e.g., Cistus spp., Genista spp.), pine plantations (Pinus pinaster), olive groves and other permanent crops, and arable cropland and pastures. Agriculture is mostly extensive, except in the Vilariça valley where there is intensive cultivation of grapes, olives, almonds, other fruits and vegetables. Throughout the region, human occupation is low (12–22 inhabitants/km2) and largely concentrated in a few larger villages.

Field sampling

Bats were captured in 2016 and 2017, aiming to obtain at least 10 individuals from all species known to occur in the region. Sampling was restricted to the period between May and July to reduce seasonal variation in prey availability, but also to match the bat's reproductive period and the activity peak of many insect pests. In total, 445 individual bats belonging to 19 different species were live trapped using mist nets placed at either the entrance of roosts, feeding perches or foraging areas. Individual bats were placed in clean cotton bags, from where guano pellets were collected. In addition, we collected 150 fresh fecal pellets from known roosts of species that could not be sampled or were undersampled by mist‐netting. Fecal pellets were stored in 2 mL tubes containing silica beads and stored at −20°C until further processing.

Laboratory procedures

Up to three fecal pellets from each bat sampled and all pellets collected in roosts were individually processed and the DNA extracted, corresponding to 1,284 individual guano pellets extractions. The number of individuals and pellets analyzed was chosen to maximize accuracy in diet estimation while controlling for costs (Mata et al. 2019). Extractions were made using a custom protocol that consisted of an initial incubation period using a lysis buffer (0.1 M Tris–HCl, 0.1 M EDTA, 0.01 M NaCl, 1% N‐lauroylsarcosine, pH 7.5–8; Maudet et al. 2002), followed by inhibitor removal using Inhibitex tablets (Qiagen, Hilden, Germany), cell lysis, DNA precipitation and washing using E.Z.N.A. Tissue Kits (Omega Bio‐Tek, Norcross, Georgia, USA). The extraction protocol was started by adding one fecal pellet to 800 µL of lysis buffer. Samples were homogenized with a spatula, vortexed, and left in a dry bath at 70°C for 30 min. Afterwards, samples were short‐spun and up to 700 µL of supernatant was transferred to a new tube containing one‐quarter of an Inhibitex tablet. Samples were then vortexed for 1 min and centrifuged at 12,000 x g for 30 s. Up to 500 µL of supernatant was transferred to a new tube and 25 µL of OB Protease was added. The remaining steps followed the kit recommendations, except that DNA was eluted two times in 50 µL into different extracts. DNA was extracted in batches of 23 samples plus one negative control in which no fecal pellet was added. Extracted DNA was distributed in 96‐well plates where the last well was left empty for as a PCR negative control.

Prey DNA was independently amplified using the ZBJ‐ArtF1c‐R2c and FwhF2‐R2n COI primer sets (Zeale et al. 2011, Vamos et al. 2017), modified to contain Illumina adaptors. ZBJ is a commonly used primer set for diet analysis of insectivorous predators that was specifically designed to amplify the most important prey of bats such as Lepidoptera and Diptera, while avoiding the amplification of the bats themselves (Zeale et al. 2011). Due to its low level of degeneracy, it is often positively biased toward these insect groups, showing a high level of mismatch to many taxa of other orders and therefore their poor amplification (e.g., Clarke et al. 2014). Fwh2 was designed to amplify freshwater invertebrates (Vamos et al. 2017), but performs quite well with terrestrial arthropods in general, being able to amplify a wide range of invertebrates, from spiders to springtails, including many non‐lepidopteran insects such as ants, beetles and true bugs (Elbrecht et al. 2019). PCR reactions consisted of 5 µL of Qiagen Multiplex Master Mix, 0.3 µL of each 10 nmol/L primer, 3.4 µL of water, and 1 µL of DNA extract. Cycling conditions consisted of a 15 min period at 95°C, 35 cycles of 30 s denaturation at 95°C, 30 s annealing at 45°C for ZBJ and at 52°C for Fwh2, and 30 s extension at 72°C, and a final extension period of 10 min at 72°C. Bat field identifications were confirmed by amplifying a small COI fragment using the bat specific primers SFF_145f‐351r (Walker et al. 2016), also adapted to be sequenced with an Illumina machine. This was done to validate the identity of cryptic species and of guano pellets collected from roosts, allowing the clear identification of all bat species occurring in Portugal (Rebelo et al. 2020), with the exception of Myotis myotis and Myotis blythii, which were therefore treated as a species complex due to well known issues of mitochondrial introgression (Afonso et al. 2017). PCR reactions and cycling conditions were similar to that of prey, except that MyTaq Mix (Bioline, Porto, Portugal) was used, and annealing was performed at 56°C. All PCR products were diluted 1:4 with water and further subjected to a second PCR reaction to incorporate 7‐bp long identification tags and Illumina P5 and P7 adaptors. PCR reactions and cycling conditions were similar to the first PCR except that KAPA HiFi HotStart ReadyMix (Rocher, KAPA Biosystems, Basel, Switzerland) was used and only eight cycles of denaturing, annealing and extension were performed, with annealing at 55°C. PCR products were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, California, USA) and subsequently quantified using a NanoDrop spectrophotometer and diluted to 15 nM. Purified and normalized PCR products were pooled per marker. These three libraries were then quantified individually using qPCR (KAPA Library Quant Kit qPCR Mix; Rocher) and diluted to 4 nmol/L. Finally, libraries were pooled by mixing 51 µL of Fwh2 library with 41 µL of ZBJ library and 8 µL of SFF library. This final library was sequenced using ˜30% of a lane of a HiSeq Rapid SBS Kit v2 (500 cycles) for a target of 38k, 25k and 5k reads/pellet for Fwh2, ZBJ and SFF amplicons, respectively. The aimed coverages were selected based on the differences in taxa amplified by the different markers. Fwh2 not only amplifies insects, but also some fungi and vertebrates such as bats, which leads to the loss of some reads to non‐dietary items. ZBJ amplifies mostly insects and does not amplify bat DNA. Finally, SFF is expected to amplify only the bat origin of the pellet as it was designed to not amplify the insects contained in the diet.

Bioinformatic analysis

Bioinformatic processing was done using Obitools (Boyer et al. 2016) and followed da Silva et al. (2019). Briefly, paired‐end reads were aligned, primer sequences removed, and reads collapsed into exact sequence variants (ESVs). ESVs that did not have a total read count in the entire dataset of 50 (approximately 0.0002% and 0.0001% read relative abundance in the entire ZBJ and Fwh2 dataset, respectively) were removed, as well as if their length was outside the expected range for the targeted taxa (202–208 bp for Fwh2 amplicons, 154–160 bp for ZBJ, and 202 bp for SFF). Finally, the ESV data were denoised to partially remove PCR and sequencing errors, using the command “obiclean” with an “r” level of 1, this way removing every “A” sequence that differed by 1 bp from sequence “B,” and had an absolute read count lower than “B” in every sample of the dataset. ESVs were then compared with online databases (BOLD and NCBI) and identified to the lowest taxonomic rank possible. Whenever an ESV matched several species, genus or families at similar identity levels, we selected the most inclusive taxonomic rank. Each ESV was also categorized as either being “diet” (i.e., most arthropods) or “not diet” (e.g., fungi, internal and external parasites). Samples that did not have at least 100 reads belonging to dietary items were considered to have failed and were discarded. Samples collected from roosts whose bat identity was not possible to assess using any of the primer pairs were also discarded. From each sample we further removed all taxa representing <1% of the total number of dietary reads of that sample to reduce eventual false interactions related to cross‐talk (Edgar 2016) and secondary predation (Deagle et al. 2019). After these procedures, the retained samples used for analysis had on average (mean ± SD) 11,719 ± 7,631 reads for ZBJ (range: 103 to 58,075) and 11,208 ± 8,036 reads for fwh2 (range: 99 to 32,764). More details on the number of reads and ESVs retained at each step of the analysis can be found in Appendix S1: Table S1. Finally, for bat individuals in which more than one pellet was analyzed, the taxa found in the different pellets were combined. The same was done for the taxa found using the two primer sets.

Analysis of bat‐pest interactions

Bat‐pest interaction patterns were described using metrics considered useful to inform conservation biocontrol (Table 1). First, we only considered prey taxa identified at the species level or that equally matched two or three species of the same genus (species complex) with high support (>99% identity). We did not consider taxa identified at higher taxonomic levels, as there would be no possible way of assessing their pest status. From the retained taxa, we identified agricultural or forest insect pests by a thorough literature survey (Appendix S2). We assumed an inclusive classification, taking as pests all species considered as such in at least a bibliographic source, even if it was not considered particularly damaging in our region. Then, we evaluated the reliability of our dataset to describe bat‐pest interactions, by building rarefaction curves based on Hill numbers for bat and pest species, and for bat‐pest interactions, using the function “iNEXT” from the package iNEXT with 1,000 bootstraps (Hsieh et al. 2016). We built both richness and sample coverage curves, considering sampling‐unit‐based incidence data (Hsieh et al. 2016).

Overall interaction patterns were described from metrics often used in molecular dietary studies (da Silva et al. 2019, Mata et al. 2019, Table 1). We computed the total number of pest species recorded in the diet of each bat species, and used the frequency of interaction (FI) to evaluate the relative abundance of pests in bat diets. We also computed the frequency of occurrence (FO) to evaluate how prevalent was pest consumption across individuals of each species. To compare the diversity of the consumed prey pests across bat species, we built rarefaction curves with 1,000 bootstraps for well sampled species (n > 10), and compared the observed value of consumed pests considering the sample size of the least sampled bat species (n = 13). For statistical significance, we considered the overlap of the 95% confidence intervals. Finally, we used a metric based on Pianka's niche overlap to compute how distinctive was pest consumption by a species in relation to all other bat species (Table 1). Overlaps were calculated with the package spa, using function “niche.overlap” with the method “pianka.” Species with less than 10 sampling units were discarded, because smaller sample sizes may produce inaccurate FO estimates (Mata et al. 2019) and therefore poor overlap estimates.

The interaction network linking each bat and pest species was built and analyzed with the R package bipartite, except otherwise indicated, and we used the number of samples at which each pest species was found to index the strength of interactions. Network maximum modularity (Q) was estimated using the function “metaComputeModules” in the package bipartite, with the “Beckett” method and 10,000 replicates (Beckett 2016). Network nestedness was calculated using the function “networklevel” with the index “weighted NODF” (Almeida‐Neto and Ulrich 2011). Network specialization level was calculated using the H 2′ index (Blüthgen et al. 2006). Statistical significance of the three network metrics was assessed by comparing the values observed with those obtained in 1,000 null models built with the function “nullmodel,” and the method “vaznull” that randomizes matrices with the same dimensions and connectivity as the initial web (Vázquez et al. 2007). As in this model interactions are assigned proportionally to a species' relative abundance, species with greater relative abundance have a higher probability of being assigned an interaction than rarely interacting species. This allowed us to test whether the modularity, nestedness and specialization of the predator‐pest interaction network was simply a consequence of the underlying species abundances or not.

The role of each bat species in the network was first assessed by calculating its centrality, considering the normalized degree, closeness and betweenness metrics (Cirtwill et al. 2018; Table 1), using the function “specieslevel.” For closeness and betweenness, we used a weighted version which accounts for the number of interactions (Opsahl et al. 2010). We then estimated the module role of each bat species, using the within‐module degree (z) and the participation coefficient (c) (Cirtwill et al. 2018, Hackett et al. 2019, Table 1), computed with the function “czvalues” on the most modular network configuration calculated previously. To define the thresholds above which a species can be considered a module hub and/or a connector, we considered the previously built 1,000 null models and used 95% quantiles as critical c‐ and z‐values, respectively (Dormann and Strauss 2014). Finally, we evaluated whether there were particular bat species specializing in pests that few other bats consumed, using the standardized specialization index, d′ (Dormann 2011), computed with function “specieslevel” (Table 1).

Finally, to explore the relative importance of different bat species to avoid pest release (i.e., pest species becoming predator‐free), we used the function “second.extinct” to build extinction curves (Memmott et al. 2004). These curves should be taken as a first and simplistic approximation, because they take the strong assumptions that before extinction the predator species actually suppressed pest species with which it interacted; that the predator needs to go fully extinct for the pest to be released; and that after extinction there is no rewiring of the network, i.e., the predators remaining in the community do not establish new interactions with the released pests. Considering these assumptions, to evaluate the relative importance of species of conservation concern, we first constructed a curve by sequentially removing bat species from highest to least conservation concern, according to their conservation status in Portugal (Cabral et al. 2005). Then, we compared this curve with that obtained from random extinctions (Memmott et al. 2004), assuming that major deviations between curves would indicate that species of conservation concern were more (or less) important than expected by chance. To estimate curves, we used 10,000 randomizations, and calculated the average number and 95% confidence interval, of predated pest species for each extinction step. In the first curve, extinctions occurred at random within conservation categories, and in the second, extinctions were completely random. To evaluate the most important species to avoid pest release, we computed a deterministic curve with extinctions occurring from species with lowest to highest degree.

Results

Overall dietary patterns

DNA metabarcoding yielded 672 interactions between 19 bat species and 132 insect pest species (Appendix S1: Table S1). Cumulative curves for richness indicated that all bat species were sampled, while only 61% of insect pest species and 41% of bat‐pest interactions were recorded (Appendix S1: Fig. S1). However, sample coverage was 91% for pests and 70% for interactions (Appendix S1: Fig. S1), indicating that few interactions with unsampled pest species were left to describe.

Bat species consumed on average 16.6 ± 10.0 (SD) pest species, with >25 species in the diets of common pipistrelle Pipistrellus pipistrellus, greater horseshoe bat Rhinolophus ferrumequinum, Schreiber's bent‐winged bat Miniopterus schreibersii, and Mediterranean horseshoe bat Rhinolophus euryale (Table 2). Although variation in the number of pest species consumed was correlated with the number of individuals analyzed per bat species (R = 0.52, P < 0.05), rarefaction analysis showed that for similar sample sizes there were some significant differences in pest species richness across bat species (Appendix S1: Fig. S2). Each pest species was consumed by 2.4 ± 2.2 (SD) bat species, with maximum values for marsh crane fly Tipula oleracea (13), pine processionary Thaumetopoea pityocampa (11), silver‐Y moth Autographa gamma (10), turnip moth Agrotis segetum (9), lesser yellow underwing moth Noctua pronuba (9), and angle shades moth Phlogophora meticulosa (8) (Appendix S1: Table S1). The FI of bats with pests varied widely (Table 2), with maximum values (>30%) for gray long‐eared bat Plecotus austriacus, European free‐tailed bat Tadarida teniotis, brown long‐eared bat Plecotus auritus, and M. schreibersii (Table 2). In most bat species, most individuals (>50%) interacted with at least one pest, with maximum values in P. austriacus. Diet distinctiveness between species was generally high (0.75 ± 0.08, 0.68–0.94) (Table 2).

Table 2

Dietary metrics, network metrics, and modularity patterns reflecting the potential contribution of individual bat species to insect pest control in a multifunctional landscape.

Species	Dietary metrics				Species roles in the network						Module
	R	FI	FO	DD	ND	BC	CC	z	c	d′
Barbastella barbastellus	9	18%	77%	0.737	0.068	0	0.036	−1.441	0	0.328	1
Miniopterus schreibersii	29	33%	86%	0.683	0.22	0.163	0.063	1.320 †	0.093	0.345	1
Plecotus auritus	15	37%	88%	0.694	0.114	0.039	0.051	−0.791	0.203	0.345	1
Plecotus austriacus	24	45%	97%	0.698	0.182	0.073	0.067	0.526	0.261	0.367	1
Rhinolophus euryale	26	25%	89%	0.73	0.197	0.112	0.06	0.513	0.322	0.415	1
Tadarida teniotis	19	39%	77%	0.682	0.144	0.051	0.061	−0.126	0.091	0.337	1
Eptesicus isabellinus	16	22%	69%	0.739	0.121	0.084	0.047	−0.322	0.457	0.396	2
Pipistrellus kuhlii	18	15%	76%	0.696	0.136	0.054	0.049	−0.055	0.524	0.382	2
Rhinolophus hipposideros	19	13%	63%	0.822	0.144	0	0.036	1.377 †	0.189	0.517	2
Myotis daubentonii	9	7%	32%	0.713	0.068	0.004	0.063	−0.548	0.374	0.424	3
Rhinolophus ferrumequinum	32	27%	79%	0.687	0.242	0.228	0.07	1.154	0.363	0.422	3
Myotis myotis/blythii	16	17%	47%	0.833	0.121	0.066	0.032	0.707	0.099	0.566	4
Myotis emarginatus	8	7%	21%	0.936	0.061	0	0.024	–	0.173	0.673	5
Myotis escalerai	16	11%	56%	0.836	0.121	0.002	0.038	−0.383	0.221	0.487	6
Pipistrellus pipistrellus	39	14%	65%	0.816	0.295	0.125	0.056	1.135	0.278	0.586	6

Bold characters indicate values within the 1st quantile. Species with <10 sampled individuals were excluded. Acronyms as in Table 1.

Species above the critical threshold to be considered a module hub.

Bat‐pest interaction network structure

The bat‐pest interaction network was modular (Q = 0.48, P‐value < 0.0001), with six modules (Fig. 1). The largest module included six bats and 48 pests, with most interactions involving M. schreibersii (22.0%), R. euryale (21.2%), P. austriacus (20.8%), and T. teniotis (17.0%), with N. pronuba (14.4%), A. segetum (12.9%), T. pityocampa (9.1%), A. gamma (8.0%), and pearly underwing moth Peridroma saucia (7.6%). The other large module included four bats and 18 pests, with most interactions involving the lesser horseshoe bat Rhinolophus hipposideros (33.9%), Kuhl's pipistrelle Pipistrellus kuhlii (33.9%), and isabelline serotine bat Eptesicus isabellinus (29.0%), with olive moth Prays oleae (21.0%), gorse shield bug Piezodorus lituratus (16.1%), and black garden ant Lasius niger (12.9%). Two other modules had three bat species each, and 29 and 31 pest species, respectively, with the first involving mainly interactions of R. ferrumequinum (60.6%) with Tipula oleracea (53.0%) and beech moth Cydia fagiglandana (12.1%), and the second interactions of P. pipistrellus (68.6%) with meadow spittlebug Philaenus spumarius (19.8%) and seedcorn maggot Delia platura (16.3%). Finally, there was one module with just two bats and 10 pests, with most interactions involving M. myotis/blythii (91.3%) with mole cricket Gryllotalpa vineae (39.1%) and the great green bush‐cricket Tettigonia viridissima (13.0%), and another module including only the Geoffroy's bat Myotis emarginatus and six pests, with most interactions involving green oak tortrix Tortrix viridana (30.0%). Network nestedness was lower than expected by chance (wNODF = 8.9, P < 0.0001), while network specialization (H 2′) was higher than expected (H 2 = 0.36, P < 0.0001), but still relatively low.

Fig. 1

Bat and pest species interaction network (a) and matrix (b). Nodes (species) and edges (interaction links) are colored according to the modules to which they were assigned. Gray edges represent interactions outside the species respective module. Node, edge and label sizes are proportional to the number of observed interactions in the network. Each bat species is represented by an abbreviation of its binomial using the first letter of the genus and the first three letters of the specific epithet (Table 2).

Bat species roles in the interaction network

Miniopterus schreibersii and R. ferrumequinum showed the highest centrality, considering their high degree, betweenness, and closeness (Table 2). P. pipistrellus and R. euryale had a high degree and betweenness, but low closeness, while P. austriacus and Daubenton's bat M. daubentonii had high closeness, but low values of the two other metrics. Regarding module‐based roles, the highest values of within‐module degree (z) were found for R. hipposideros, M. schreibersii, R. ferrumequinum, and P. pipistrellus, each belonging to a different module. However, only the former two species were above the threshold (z > 1.28) to be considered module hubs (Table 2). The participation coefficients were all well below the critical threshold for a species to be considered a connector (c > 0.65), which together with the relatively low values of within‐module degree categorized all but the two module hubs as peripheral species (Table 2). The average specialization index d′ was relatively low (mean ± SD: 0.46 ± 0.11), with the highest values (>0.50) for M. emarginatus, P. pipistrellus, M. myotis/blythii and R. hipposideros (Table 2).

Vulnerability to extinction

The extinction curves suggested that species of conservation concern were more important to avoid pest species release than expected by chance (Fig. 2). The deterministic extinction curve revealed that the three species with highest degree (P. pipistrellus, R. ferrumequinum and R. euryale) were sufficient to ensure that ˜60% of pest species had at least one bat predator. This proportion raised to ˜75% by adding the next three species with highest degree (M. schreibersii, P. austriacus, and T. teniotis). From these six species, two are vulnerable (R. ferrumequinum and M. schreibersii), one is critically endangered (R. euryale) and one is data deficient (T. teniotis) in Portugal.

Fig. 2

Effect of bat species extinction on the percentage of pest species with at least one bat predator. Curves were built assuming three different scenarios based on the number of pest interactions (abundance), randomly (random), and by decreasing level of conservation status (conservation status).

Discussion

We provide a framework combining DNA metabarcoding and ecological network analysis to identify small vertebrate species potentially providing key services of pest control in complex multifunctional landscapes. The application of this framework to a Mediterranean landscape revealed with unprecedented detail a highly complex trophic network involving 19 bat species feeding on over 130 insect pests, with some generalists and a few specialist predators consuming most pests. We suggest that this approach can be widely used to inform conservation biocontrol, thereby contributing to enhance the potential for multifunctional landscapes integrating biodiversity conservation with the sustainable production of food and fiber.

Bat dietary patterns in multifunctional landscapes

Our study revealed the complexity of trophic interactions occurring in multifunctional landscapes between highly mobile vertebrate predators and a wide range of insect pests. The number of pest species recorded in bat diets was very high (132), with rarefaction curves suggesting that the true number may still be ˜40% higher. These results agree with previous research showing that individual bat species can forage on a wide range of pests (Aizpurua et al. 2018, Cohen et al. 2020), and that individual pest species may be preyed upon by several bat predators (Garin et al. 2019). However, the diversity of pests was much higher in ours than in previous studies, possibly because we sampled a wide range of bat species with different habitat preferences and foraging strategies, which in turn may contribute to increase the range of prey consumed by the predator community. Moreover, bats in our area probably had access to a wide range of agricultural and forest pest species, due to the coexistence at the landscape scale of multiple land uses. The high number of pest species recorded was possibly also due to our large sampling effort compared with previous studies (e.g., Aizpurua et al. 2018, Cohen et al. 2020), with over 1,200 guano pellets individually analyzed, and the use of two different primer sets, one of which (Fwh2) had a high degeneracy level and could therefore capture a wide diversity of prey. Whatever the reasons, these results pointed out the need for community‐level studies to understand predator‐pest interactions, because focusing on a limited number of predators and pests may miss the full array of potentially relevant interactions occurring at the landscape scale.

Although all bat species in our study preyed on pests, the number of species consumed ranged from <10 in M. emarginatus, B. barbastellus and M. daubentonii, to >25 in R. euryale, M. schreibersii, R. ferrumequinum, and P. pipistrellus. Although this result was influenced by variation across bat species in the number of individuals analyzed, there were still significant differences after controlling for sample size effects. This means that at least part of the variation was probably real, and may be a consequence of some bat species feeding preferentially on highly diverse pest groups such as moths, having generalist feeding strategies or foraging in habitats that favor the consumption of a wide range of prey species. There was also wide interspecific variation in the FI with pest species, with particularly high values (>30%) for M. schreibersii, P. auritus, T. teniotis and P. austriacus. Although FI does not translate directly into food intake, it provides an indication that pest species were far more consumed by some bat species than others. This probably reflect a matching between prey availability and dietary preferences, with higher FIs for bat species feeding on abundant pests (Aizpurua et al. 2018, Costa et al. 2020). These results highlight major differences across bat species in their interaction with pests, which may affect their role in biocontrol.

We also observed wide variation in the consumption of different prey, with ˜8% of the pest species being responsible for nearly 50% of all bat‐pest interactions. All these pests were eaten by at least five bat species, with up to 13 species feeding on T. oleracea, 11 on T. pityocampa, and 10 on A. gamma. Some of the heavily eaten pests are known to cause significant economic damages, including Prays oleae on olive orchards, T. pityocampa on pine plantations, and A. segetum and P. saucia on a variety of crops (Aizpurua et al. 2018). P. spumarius is a vector of the plant pathogenic bacterium Xylella fastidiosa, which can have strongly negative economic impacts (Cornara et al. 2017). Heavy predation on these pests is probably due to their high availability across the landscape, although this would need to be confirmed with data on species abundances. Moreover, these are probably pests fitting in the trophic niche of many generalist and specialist bat species (e.g., Mata et al. 2016, Garin et al. 2019, Cohen et al. 2020), thereby favoring their widespread consumption. Overall, these results support the view that bat species from different trophic guilds can concentrate feeding on abundant pest species, which can represent valuable food resources (McCracken et al. 2012, Aizpurua et al. 2018, Garin et al. 2019).

Learning from predator‐pest networks to inform conservation biocontrol

The analysis of predator‐pest interactions provided useful information to identify bat species potentially playing key roles in pest regulation in our landscape (Table 1), and that might therefore be priority targets for conservation biocontrol (Shields et al. 2019). Dietary data showed that R. euryale, R. ferrumequinum, M. schreibersii, P. pipistrellus, T. teniotis, P. auritus and P. austriacus may be particularly important because they feed on a wide range of pest species, and/or incorporate a large proportion of pests in their diets (Aizpurua et al. 2018, Garin et al. 2019, Cohen et al. 2020). These data also showed high diet distinctiveness, indicating that a wide range of bat species may be needed to regulate the full complement of pest species. However, more detailed information could be obtained by looking at the interaction networks in which these species were integrated.

The organization of our network in six modules indicated that there were sets of bat and pest species more likely to interact with each other (Dormann and Strauss 2014). This implies that at least some bat species representative of each module need to be targeted by conservation biocontrol to control most pest species, and to contribute for network stability by maintaining its modularity (Fortuna et al. 2010). Modularity information can also help fine tuning conservation biocontrol of specific pests that are regionally relevant. For instance, considering potentially damaging and heavily predated pests in our landscape, conservation biocontrol of A. segetum, T. pityocampa and P. saucia should probably target major pest predators such as M. schreibersii, P. austriacus and T. teniotis, while that of P. oleae should target R. hipposideros and P. kuhlii. Similarly, conservation biocontrol of P. spumarius should target bat species from a small module that includes P. pipistrellus. This pattern of modularity, together with a high connectivity between bats and pests, was probably responsible for the lack of nestedness in our network (Fortuna et al. 2010), which implies that eventual biocontrol services provided by specialists were not redundant to those provided by generalists. Modularity was probably also responsible for network specialization (H 2′) to be higher than expected by chance, although its relatively small value indicates that many bat species were connected to several pest species.

Some species appeared to have keystone roles in the interaction network (González et al. 2010, Cirtwill et al. 2018), and may therefore be key targets for conservation biocontrol. These included species with high degree and betweenness centrality such as M. schreibersii, R. ferrumequinum P. pipistrellus, and R. euryale, which may be essential to ensure coherence in the network, and therefore to maintain its structure and stability (González et al. 2010). The former two species, plus P. austriacus and M. daubentonii, also had the highest closeness centrality and are therefore were important to the overall connectivity of the network (González et al. 2010). Based on values of within‐module degree, there also appeared to be one keystone species in each of the four largest modules, including R. hipposideros, M. schreibersii, R. ferrumequinum and P. pipistrellus, although only the former two species were considered module hubs. This suggests that each of these species may play a unique role at regulating particular groups of pest species. The specialization index further highlighted species with unique roles in pest consumption, including M. emarginatus, P. pipistrellus, M. myotis/blythii, and R. hipposideros. Finally, extinction curves showed that just six bat species were sufficient to ensure that ˜75% of pest species had at least one bat predator. Three of these species are of conservation concern in Portugal, highlighting the potential importance of endangered species to pest control.

Combining the different metrics analyzed, our results suggested that M. schreibersii, R. hipposideros, R. ferrumequinum and P. pipistrellus may be critical targets for conservation biocontrol in our landscape, as they feed heavily on pests, together they cover most pest species, and play important roles at the network and module levels. To these should be added some specialized species such as M. myotis/blythii and M. emarginatus that feed on pests rarely consumed by other bats. Finally, additional species are important because they also feed heavily on many pest species and provide redundancy to pest control services, including R. euryale, T. teniotis and P. austriacus. Despite this identification of priority species for conservation biocontrol, the overarching objective of maintaining the integrity of the overall bat community should not be neglected, not only for assuring the persistence and stability of pest control services, but also to fulfill the objective of conserving biodiversity in multifunctional landscapes.

Limitations and future developments

The application of our framework to inform conservation biocontrol needs to account for some limitations and caveats, part of which are inherent to DNA metabarcoding. One potential problem is secondary ingestion, as metabarcoding may detect DNA from species represented in predator gut contents (e.g., da Silva et al. 2019, Deagle et al. 2019). As such, we may have recorded some pest species consumed by arthropod predators (e.g., arachnids) rather than directly by bats. In general, however, this problem is unlikely to be serious, because DNA derived from gut contents tends to be more degraded and to occur in smaller quantities than that from prey ingested directly, and so the latter are far more likely to be recovered in larger quantities during PCR amplification. Still, we recommend the approach taken in our study to filter out species with a very low proportion of reads (<1%) in each sample, thereby minimizing the detection of secondary ingestion. Moreover, we recommend that metabarcoding results should be checked against ancillary ecological data, to check for their consistency with previous information on the diet and foraging behavior of studied predators (da Silva et al. 2019, Deagle et al. 2019), as observed in our study. Conversely, however, secondary ingestion may provide novel opportunities to further understand the network of trophic interactions between bats and arthropod pests, as plant DNA possibly originated from the guts of prey can be recovered in vertebrate fecal samples (da Silva et al. 2019), which in turn might be used to obtain information on economically important crops that pest species were feeding upon.

Metabarcoding is also unable to provide reliable information on the number, or biomass, of each prey consumed, making it difficult to infer interaction strength between predators and prey. Our framework circumvented this problem by considering the FI of each bat and pest species, which indicates how often a predator consumes a prey and can be used as a surrogate of interaction strength. This approach provides more complete information than just considering interaction presence/absence, but it requires large samples sizes to accurately estimate frequencies of interaction (Mata et al. 2019). Therefore, the application of our framework should be based on a comprehensive sampling effort, which in our case involved large numbers of individuals per bat species and the analysis of multiple pellets per individual. Still, this sampling effort would probably have to be increased to provide more accurate estimates for rarer predators and pest species ingested less frequently (Mata et al. 2019). However, this problem should have a smaller effect on estimates for common predators and pests, which are likely to be the most relevant for conservation biocontrol.

The application of our framework may also be limited by the sampling effort needed to produce detailed ecological networks, because several hundred samples are usually required to capture the diversity of communities, and even more to fully describe species interactions in those communities (Chacoff et al. 2012, Jordano 2016). Therefore, it is fundamental to estimate sampling completeness, which in our case showed that the diversity of pests and bat‐pest interactions were not fully recovered. Yet, sample coverage suggested that we mainly missed rare pest species and interactions, which are probably less relevant for conservation biocontrol. Moreover, the effects of undersampling on network metrics such as modularity, nestedness and specialization seem to be limited when at least 30% of the species are recorded (Blüthgen et al. 2006, Rivera‐Hutinel et al. 2012, Costa et al. 2016), which suggests that robust inferences on network structure can be made even when sampling completeness is lower than in our study. Regarding species roles in networks, the impacts of undersampling are much less understood, although they need to be considered because metrics such as centrality may be sensitive to the number of interactions detected, which in turn is influenced by sample size (Fründ et al. 2016). To minimize this problem, in our study we did not consider species with low sample sizes (≤10), while for the others we found that correlations between the number of individuals analyzed and species roles' metrics were always small (all R < 0.36). Overall, therefore, we suggest that meaningful information for conservation biocontrol may be obtained from ecological networks based on sampling efforts comparable with ours, although increased sample sizes may be needed in applications requiring information on rarer predator species and/or involving rare interactions with pest species.

Finally, although our approach can be used to identify key candidate species potentially involved in pest control, it is not sufficient to quantify the actual delivery of such service. For this, far more detailed information would be needed, including that on the actual predation pressure on pest species, which in turn requires data on predator population sizes and on the number or biomass of each pest consumed by each predator (Feit et al. 2019). Moreover, attention should be given to potential disservices associated with intraguild predation, because top predators feeding on arthropod predators of agricultural and forest pests may reduce net predation pressure on pests (e.g., Garcia et al. 2020, Garfinkel et al. 2020, Olimpi et al. 2020). Addressing this issue would require for instance the analysis of comprehensive networks including also beneficial predators, based on dietary data for both top and mesopredators. In any case, fully demonstrating that predation has an actual effect requires experimental approaches, whereby changes in pest population sizes or damages are evaluated against artificial manipulation of predation pressure (Maine and Boyles 2015, Maas et al. 2016). These approaches provide the gold standard to show the actual importance of a predator to pest control, but we argue that they may be generally unfeasible in landscapes such as ours, where there are many crops and forest production systems, and therefore multiple predators interacting with multiple pests across habitats. In these circumstances, a relatively simple approach such as ours represents an important advance toward conservation biocontrol, providing a knowledge basis to produce adaptive management prescriptions that can later be refined once more detailed data become available.

Conclusions

This study shows the power of combining DNA metabarcoding and ecological network analysis to guide ecosystem management, by providing a detailed understanding of interactions between entire communities. This possibility has already been raised (e.g., Evans et al. 2016, Bohan et al. 2017, Clare et al. 2019), but few studies, if any, have actually used this approach to address practical management problems. Here we provide such practical application, showing how the approach can help to identify bat species that may have key roles for pest control in complex multifunctional landscapes. Although this approach cannot demonstrate actual pest control by these species, and the quality of inferences are limited by the quality and quantity of data used to build the network, we are convinced that it can bring useful information to guide conservation biocontrol that would be difficult to gather otherwise. For instance, our study strongly suggests that some of the bats potentially critical for biocontrol are also endangered species, and so conservation efforts targeted at these species may have wider benefits for pest control at the landscape scale. We believe that approaches such as ours may find wide applicability in the near future, as costs of metabarcoding are ever decreasing, and technical problems are being solved (e.g., da Silva et al. 2019, Deagle et al. 2019, Mata et al. 2019, Piñol et al. 2019). This in turn may help to improve the management of multifunctional landscapes, where biodiversity conservation and the delivery of ecosystem services are critically dependent on the interactions established between multiple species across trophic levels.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.