A molecular analysis of the diet and biocontrol potential of spiders in cereal crops - Dataset

Methodology

Fieldwork

Money spiders (Araneae: Linyphiidae) and wolf spiders (Araneae: Lycosidae) were visually located along transects in two adjacent barley fields at Burdons Farm, Wenvoe in South Wales (51°26'24.8"N, 3°16'17.9"W) and collected from occupied webs and the ground, between April and September 2018. Each belt transect was adjacent to a randomly selected crop tramline and they were distributed across the entire field and ran its length. The areas searched were 4 m2 quadrats at least 10 m apart and all observed linyphiids and lycosids were collected. Spiders were taken from 64 locations along the aforementioned transects. Spiders were placed in 100 % ethanol using an aspirator, regularly changing meshing to limit potential cross-contamination between spiders introduced by alternative techniques such as pitfall or suction sampling1,2. Linyphiids occupying webs were prioritised for collection, but ground-hunting linyphiid spiders were collected when occupied webs were scarce. For each individual money spider that was taken from its web, the height of the web from the ground and its approximate dimensions were recorded, the latter calculated as approximate web area. Spiders were taken to Cardiff University, transferred to fresh ethanol, adults identified to species-level and juveniles to genus, and stored at -80 °C in 100 % ethanol until subsequent DNA extraction. To obtain data on local prey density, 4 m2 of ground and crop stems was suction sampled for approximately 30 seconds at each quadrat from which spiders were collected, with the collected material emptied into a bag and any organisms immediately killed with ethyl-acetate. Suction sampling used a ‘G-vac’ modified garden leaf-blower, known to have a greater air velocity and capture more taxa relevant to this study, such as spiders, hymenopterans and thrips, when compared to commonly-used Vortis™ samplers3. All vacuum-sampled material was frozen for storage before sorting into 70 % ethanol in the lab.

All invertebrates were identified to family level due to the restriction of many of the metabarcoding-derived dietary data to this level, and the difficulty associated with finer taxonomic resolution of many damaged specimens. Exceptions included springtails of the superfamily Sminthuroidea (Sminthuridae and Bourletiellidae, which were often indistinguishable following suction sampling and preservation due to the fine features necessary to distinguish them) which were left at super-family, mites (many of which were immature or in poor condition) which were identified to order level and wasps of the superfamily Ichneumonoidea (which were identified no further due to obscurity of wing venation due to damage). Identifications were carried out under an Olympus SZX7 stereomicroscope using morphological keys4–10.

Extraction and high-throughput sequencing of spider gut DNA

Given their prevalence in field collections, dietary analysis was carried out for the linyphiid genera Erigone, Tenuiphantes, Bathyphantes and Microlinyphia (Araneae: Linyphiidae), and the Lycosidae genus Pardosa. Spiders were transferred to and washed in fresh 100 % ethanol to reduce external contaminants prior to identification via morphological key5. Abdomens were removed from spiders and again transferred to and washed in fresh 100 % ethanol. DNA was extracted from the abdomens via Qiagen TissueLyser II and DNeasy Blood & Tissue Kit (Qiagen) as per the manufacturer protocol, but with an extended lysis time of 12 hours to account for the complex and branched gut system in spider abdomens11.

For amplification of DNA, two primer pairs were used. BerenF-LuthienR12 amplified a broad range of invertebrates including spiders, and TelperionF-LaureR, amplified a range of invertebrates but fewer spiders (modified from TelperionF-LaurelinR12 via one base-pair change; 5’-ggrtawacwgttcawccagt-3’). These two primer pairs amplified 314 bp (BerenF-LuthienR) and 302 bp (TelperionF-LaureR) regions of the COI gene. Primers were labelled with unique 10 bp molecular identifier tags (MID-tags) so that each individual had a unique pairing of forward and reverse for identification of each spider post-sequencing. PCR reactions of 25 µl volumes contained 12.5 µl Qiagen PCR Multiplex kit, 0.2 µmol (2.5 µl of 2 µM) of each primer and 5 µl template DNA. Reactions were carried out in the same thermocycler, optimised via temperature gradient, with an initial 15 minutes at 95 °C, 35 cycles of 95 °C for 30 seconds, the primer-specific annealing temperature for 90 seconds and 72 °C for 90 seconds, respectively, followed by a final extension at 72 °C for 10 minutes. BerenF-LuthienR and TelperionF-LaureR used annealing temperatures of 52 °C and 42 °C, respectively.

Within each PCR 96-well plate, 12 negative controls (extraction and PCR), 2 blank controls and 2 positive controls were included (i.e. 80 samples per plate), based on Taberlet et al.13. Positive controls were mixtures of invertebrate DNA comprised of non-native Asiatic species in four different proportions (Supplementary Table 1) and blanks were empty wells within each plate to identify tag-jumping into unused MID-tag combinations. PCR negative controls were DNase-free water treated identically to DNA samples. A negative control was present for each MID-tag to identify any contamination of primers. All PCR products were visualised in a 2 % agarose gel with SYBRSafe (Thermo Fisher Scientific, Paisley, UK) and placed in categories based on their relative brightness. The concentration of these brightness categories was quantified via Qubit dsDNA High-sensitivity Assay Kits (Thermo Fisher Scientific, Waltham, MA, USA) with at least three representatives of each category per plate. The PCR products were then proportionally pooled according to these concentrations. Each pool was cleaned via SPRIselect beads (Beckman Coulter, Brea, USA), with a left-side size selection using a 1:1 ratio (retaining ~300-1000 bp fragments). The concentration of the pooled DNA was then determined via Qubit dsDNA High-sensitivity Assay Kits and pooled together into one library per primer pair. Library preparation for Illumina sequencing was carried out on the cleaned libraries via NEXTflex Rapid DNA-Seq Kit (Bioo Scientific, Austin, USA) and samples were sequenced on an Illumina MiSeq via a V3 chip with 300-bp paired-end reads (expected capacity ≤25,000,000 reads).

Bioinformatic analysis

The Illumina run generated 11,165,405 and 10,959,010 reads for BerenF-LuthienR and TelperionF-LaureR, respectively, which were quality-checked and paired via FastP14  to retain only sequences of at least 200 bp with a quality threshold of 33, resulting in 10,561,874 and 9,355,112 paired reads. The paired reads were demultiplexed and assigned to their respective spider sample according to their MID-tags via the “trim.seqs” command in Mothur v1.39.515, leaving 7,854,610 and 7,437,929 reads with exact matches to the primer and MID-tags.

Replicates were removed, and denoising and clustering to zero-radius operational taxonomic units (ZOTUs; clustered without % identity to avoid multiple species represented within a single operational taxonomic unit (OTU)) completed via Unoise3 in Usearch1116. The resultant sequences were assigned a taxonomic identity from GenBank via BLASTn v2.7.117 using a 97 % identity threshold18. The BLAST output was analysed in MEGAN v6.15.219. Where the top BLAST hit, determined by lowest e-value, was resolved at a higher taxonomic level than species-level, the results were checked; where possibly erroneous entries were preventing species-level assignment (e.g. poorly-resolved identifications on GenBank), finer resolution was assigned based on the next-closest match. Where ZOTUs were assigned the same taxon, these were aggregated.

Data clean-up followed the protocol described by Drake et al.20. The maximum value for a ZOTU present in blank or negative controls was identified and subtracted from all read counts for that ZOTU to remove background contaminants. Simultaneously, known lab contaminants (e.g. German cockroach Blattella germanica), artefacts and errors of the sequencing process, unexpected reads in positive controls and positive control taxon reads in dietary samples were identified. These were calculated as a percentage of their respective sample’s read count and any read counts lower than the highest of these percentages for their respective sample were removed to eliminate additional instances of contamination. These thresholds were defined as 0.38 % and 0.39 % for BerenF-LuthienR and TelperionF-LaureR, respectively. The data from the two libraries (i.e. from each primer pair) were then aggregated together by sample and aggregated again by taxon. Non-target taxa (e.g. fungi) and instances in which predator DNA was amplified (i.e. ZOTUs with high read counts matching the individual’s morphological identity) were removed. All remaining read counts were converted to presence-absence.

Statistical analysis

All analyses were conducted in R v4.0.021. Prey species that occurred only once across all of the dietary samples were removed before further analyses to prevent outliers skewing the results, which is particularly problematic for non-metric multidimensional scaling. Spider diets were compared between variables using multivariate generalized linear models (MGLMs) via ‘manyglm’ in the ‘mvabund’ package22 with a binomial error family and Monte Carlo resampling. Model independent variables included spider genus, spider life stage (juvenile or adult, the latter defined by fully developed genitalia), spider sex and all two-way interactions between these variables.

Coarse dietary differences were visualised by non-metric multidimensional scaling (NMDS) via metaMDS in the ‘vegan’ package23 with Jaccard distance in two dimensions and 999 tries. For NMDS, outliers (usually samples containing rare taxa) were identified by plotting and subsequently removed to facilitate separation of samples and achieve minimum stress. For visualisation of the effect of categorical variables against the dietary NMDS, spider plots were created using ‘ordispider’ with ‘ggplot’ and the ‘RColorBrewer’ ‘Accent’ colour palette24. Spider diet was compared against web characteristics for spiders for which both data were available using the MGLM process outlined above, but with starting models containing only web height, web area and an interaction between the two, with the same binomial error family, but with a ‘cloglog’ link function.

All prey taxa were classified as agricultural pests, natural enemies or excluded from subsequent analyses of intraguild predation and biocontrol (Supplementary Table 2). Intraguild predation and biocontrol variables were created by counting the number of natural enemy taxa, and, separately, of agriculturally relevant “pest” taxa (taxa containing species that commonly detract from agricultural productivity; Supplementary Table 2) in each spider’s diet. These were analysed against spider genus, maturity and sex via GLMs. “Site” (denoting the 4 m2 area from which spiders were collected within fields) was initially included as a random effect in generalized linear mixed-models, but no significant effect was observed when comparing this model against a standard GLM via a likelihood ratio test of nested models using the ‘lrtest’ command in the ‘lmtest’ package25. Standard GLMs were thus used to avoid issues relating to singularity in the mixed models. The assumptions for the resultant Poisson error family GLMs were tested using the “testResiduals” function of the ‘DHARMa’ package26. Intraguild predation and biocontrol differences between significant terms were visualised using violin plots with the quartiles, median and 95 % upper limit annotated using the ‘geom_violin’ function in ‘ggplot2’.

In situ spider prey choice was analysed using network-based null models in the ‘econullnetr’ package27 with the ‘generate_null_net’ command, visually represented with the ‘plot_preferences’ command. Suction sample data were used to represent prey availability. Prior to the taxonomic prey choice analysis, an hemipteran identified no further than order level through dietary analysis was removed due to the inability to pair it to any present prey taxa with certainty. Standardised effect sizes (SES) were extracted for all comparisons for each individual spider and compared between genera, life stages and sexes using permutational multivariate analysis of variance (PerMANOVA) using the ‘adonis’ function of the ’vegan’ package with 9999 permutations and a Euclidean distance matrix to determine overall differences in prey choice. Similarity percentages (SimPer) analysis was used to identify any differences in spider group preference for specific prey taxa using the ‘simper’ function of the ‘vegan’ package with 9999 permutations.

References

1.        King, R. A. et al. Suction sampling as a significant source of error in molecular analysis of predator diets. Bull. Entomol. Res. 102, 261–266 (2012).

2.        Athey, K. J., Chapman, E. G. & Harwood, J. D. A tale of two fluids: does storing specimens together in liquid preservative cause DNA cross-contamination in molecular gut-content studies? Entomol. Exp. Appl. 6, 338–343 (2017).

3.        Zentane, E., Quenu, H. & Graham, R. I. Suction samplers for grassland invertebrates: comparison of numbers caught using Vortis and G-vac devices. Insect Conserv. Divers. 9, 470–474 (2016).

4.        Goulet, H. & Huber, J. T. Hymenoptera of the World: An Indentification Guide to Families. (Agriculture Canada, 1993).

5.        Roberts, M. J. The Spiders of Great Britain and Ireland (Compact Edition). (Harley Books, 1993).

6.        Unwin, D. A Key to the Families of British Bugs (Insecta, Hemiptera). (FSC Publications, 2001).

7.        Ball, S. Introduction to the Families of British Diptera Part 2: Key to Families. (Dipterists Forum, 2008).

8.        Barber, A. D. Key to the Identification of British Centipedes. (FSC Publications, 2008).

9.        Duff, A. G. Beetles of Britain and Ireland Volume 1: Sphaeriusidae to Silphidae. (A.G. Duff Publishing, 2012).

10.      Dallimore, T. & Shaw, P. Illustrated Key to the Families of British Springtails (Collembola). (FSC Publications, 2013).

11.      Krehenwinkel, H., Kennedy, S., Pekár, S. & Gillespie, R. G. A cost-efficient and simple protocol to enrich prey DNA from extractions of predatory arthropods for large-scale gut content analysis by Illumina sequencing. Methods Ecol. Evol. 8, 126–134 (2017).

12.      Cuff, J. P. et al. Money spider dietary choice in pre- and post-harvest cereal crops using metabarcoding. Ecol. Entomol. 46, 249–261 (2021).

13.      Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. Environmental DNA. (Oxford University Press, 2018).

14.      Chen, S., Zhou, Y., Chen, Y. & Gu, J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890 (2018).

15.      Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).

16.      Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461 (2010).

17.      Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics 10, 1–9 (2009).

18.      Alberdi, A., Aizpurua, O., Gilbert, M. T. P. & Bohmann, K. Scrutinizing key steps for reliable metabarcoding of environmental samples. Methods Ecol. Evol. 9, 1–14 (2017).

19.      Huson, D. H. et al. MEGAN Community Edition - interactive exploration and analysis of large-scale microbiome sequencing data. PLoS Comput. Biol. 12, 1–12 (2016).

20.      Drake, L. E. et al. Post-bioinformatic methods to identify and reduce the prevalence of artefacts in metabarcoding data. Authorea April 13 (2021) doi:https://doi.org/10.22541/au.161830201.18684167/v1.

21.      R Core Team. R: A language and environment for statistical computing. (2020).

22.      Wang, Y., Naumann, U., Wright, S. T. & Warton, D. I. mvabund – an R package for model-based analysis of multivariate abundance data. Methods Ecol. Evol. 3, 471–474 (2012).

23.      Oksanen, J. et al. vegan: Community Ecology Package. (2016).

24.      Neuwirth, E. RColorBrewer: ColorBrewer palettes. (2014).

25.      Zeileis, A. & Hothorn, T. Diagnostic checking in regression relationships. R News 2, 7–10 (2002).

26.      Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. (2020).

27.      Vaughan, I. P. et al. econullnetr: an r package using null models to analyse the structure of ecological networks and identify resource selection. Methods Ecol. Evol. 9, 728–733 (2018).