Evidence for nutrient-specific foraging of predators under field conditions

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 were distributed across the entire field and ran its length. The areas searched were 4 m² quadrats at least 10 m apart and all observed linyphiids and lycosids were collected. Spiders were taken from 64 randomly selected locations along the aforementioned transects. Following collection of spiders, 4 m² of ground and crop stems was suction sampled for approximately 30 seconds, 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. All material was later frozen at -20 ºC for storage before sorting in the lab. These invertebrates were collected for background population densities and not for any molecular work.

All invertebrates were identified to family level. Further identifications were not carried out due to the inability to identify some of the invertebrate groups further via the associated metabarcoding-derived dietary data (e.g., Sciaridae), and the difficulty associated with finer taxonomic resolution of many damaged or immature specimens. The only taxa not identified to family level were springtails of the superfamily Sminthuroidea (Sminthuridae and Bourletiellidae, which were often indistinguishable following suction sampling and preservation due to the fine features necessary to differentiate them) which were left at super-family, mites (many of which were immature or in poor condition, or lacked appropriate taxonomic keys) 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); in these cases, these taxonomic assignments were pooled to family-level for later analyses.

Extraction and high-throughput sequencing of spider gut content DNA

Given their prevalence in field collections, dietary analysis was carried out for the linyphiid spider genera Erigone, Tenuiphantes, Bathyphantes and Microlinyphia (Araneae: Linyphiidae), and Pardosa (Araneae: Lycosidae). Spiders were transferred to and washed in fresh 100% ethanol to reduce external contaminants prior to identification via morphological keys(1). 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 abdomens(2).

For amplification of DNA, two primer pairs were used. BerenF-LuthienR(3) amplified a broad range of invertebrates including spiders, and TelperionF-LaureR, amplified a range of invertebrates with the exception of some spiders (modified from TelperionF-LaurelinR(3) (via one base-pair change to decrease host DNA amplification; 5’-ggrtawacwgttcawccagt-3’). These two primer pairs amplified 314 bp (BerenF-LuthienR) and 302 bp (TelperionF-LaureR) regions of COI. 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, optimized 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.(4). Positive controls were mixtures of invertebrate DNA comprised of non-native Asiatic species in four different proportions (Table S1) 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 visualized in a 2 % agarose gel with SYBR®Safe (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 FastP(5)  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.5(6), 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 amplicon sequence variants (ASVs; clustered without % identity to avoid multiple species represented within a single operational taxonomic unit (OTU)) completed via Unoise3 in Usearch11(7). The resultant sequences were assigned a taxonomic identity from GenBank via BLASTn v2.7.1(8) using a 97% identity threshold(9). The BLAST output was analyzed in MEGAN v6.15.2(10). 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 ASVs were assigned the same taxon, these were aggregated.

Data clean-up followed the protocol described as optimal by Drake et al.(11). The maximum value for an ASV present in blank or negative controls was identified and subtracted from all read counts for that ASV 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. ASVs with high read counts matching the individual’s morphological identity) were removed. All remaining read counts were converted to presence-absence.

Macronutrient determination

Specimens were taken for macronutrient analysis from the same suction samples collected for invertebrate community identification. Representatives were taken from each family found in the community samples for which specimens were intact, in visually good condition and relatively clean of soil and other contaminants. If specimens were from a relatively uncommon family but unclean, soil and other surface contaminants were physically removed, and the specimen then momentarily dipped in water to remove remaining surface contaminants without greatly dislodging surface lipids. Macronutrient contents were determined following the MEDI protocol(12, 13) with minor alterations to account for the small size of most of the invertebrates processed(14). During extraction, half volumes (i.e. 500 µl) of solvents were used. For the lipid assays, 15 µl of sulfuric acid was added for a 15 min incubation, followed by only 200 µl of vanillin reagent to increase the concentration and development of analyte for more accurate readings from smaller invertebrates. Lipid and protein standard series were diluted to 50% of the concentration specified in the original protocol (i.e. 0-1 mg ml^-1). Carbohydrate assays used 140 µl of reagent with 30 min incubation at 92 °C followed by a further 30 min at room temperature. Carbohydrate standard series were diluted to 1% of the concentrations specified in the original protocol (i.e. 0-0.02 mg ml^-1) to ensure signals overcame the higher limit of detection relative to typical invertebrate carbohydrate content.

Statistical analysis

All analyses were conducted in R v.4.0.3(15). In situ spider prey choice was analyzed using network-based null models in econullnetr(16) with the ‘generate_null_net’ command. A bespoke set of functions was used alongside econullnetr to randomly generate an “expected diet” for each individual spider based on local prey communities determined via suction sampling. Macronutrient data were allocated to each dietary taxon and the mean macronutrient proportions calculated. The mean macronutrient contents were compared between expected and observed diets using a multivariate linear model (MLM) via mvabund(17) and significant differences visually represented through ternary plots using ggtern(18) and ggplot2(19). The observed mean nutrient proportions of spider diets were compared between spider genera, life stages and sexes using a MLM. To ascertain how prey choice factors into these dietary differences, the difference in macronutrient proportions between expected and observed spider diets were also compared between spider genera, life stages and sexes in a MLM.

To group taxa into tropho-species, mean macronutrient values for each taxon were first determined to prevent splitting of taxa across clusters; these were represented at the family, order and class levels to allow tropho-species assignment for families for which macronutrient content was not determined, but was at a higher level. Macronutrient values were scaled by subtracting the mean of each column from each contained value and dividing it by the column standard deviation using the ‘scale’ function. A Euclidean distance matrix was calculated using the ‘dist’ function. Hierarchical clustering of scaled macronutrient distance matrix used the ‘hclust’ function. Optimal clustering solutions were determined by comparison of Dunn’s index between methods and k values; this was calculated using the ‘dunn’ function in the “clValid” package(20) for each cluster k value above five until the Dunn index decreased, the first instance of the value preceding the decrease deemed the maximum value, thus optimal solution. Clustering solutions based on ‘average’, ‘complete’, ‘single’, ‘median’, ‘centroid’ and ‘mcquitty’ linkages were compared, and the “complete” method selected for subsequent analysis as it resulted in the smallest number of clusters (20; thus, the most efficient simplification of the taxa analyzed). Three uncommon families (present in small numbers in one community sample each, but no dietary detections) were removed from further tropho-species analyses due to the lack of class-level macronutrient data (Arionidae, Lithobiidae and Polydesmidae).

To name the tropho-species, a second clustering stage was used in which the tropho-species were grouped according to their mean macronutrient content for each of the three nutrients separately. ‘Single’ linkage clustering was found to be the optimal method for this step and created ten, seven and six groups for carbohydrate, lipid and protein, respectively. These clusters were labelled from one to the total number of clusters for each macronutrient to represent low-to-high content of that nutrient relative to other tropho-species. Names used the structure ‘CxLyPz’ to denote the relative content of each tropho-species (x, y and z replaced with the cluster number for carbohydrate, lipid and protein, respectively).

Clusters were henceforth termed ‘tropho-species’, with all taxa within a single cluster representing a single aggregated tropho-species. Heatmap dendrograms were produced using the ‘heatmap.2’ function in the ‘gplots’ package(21), with cluster colors assigned with the ‘Accent’ palette of ‘RColorBrewer’(22) and relative macronutrient content color scaling produced using the ‘viridis’ package(23). Ternary plots were produced to visualize the macronutrient content of taxa within each cluster, and differences in mean macronutrient contents between tropho-species.

Tropho-species were assigned to each taxon present in dietary and prey community samples. Where family-level macronutrient data were not obtained (usually low abundance and poor condition invertebrates or families identified in the diet that were not subsequently observed in community samples), order-level tropho-species assignment was used, or class where order-level data were not available (12 and 2 instances of uncommon taxa, respectively).

In situ spider prey choice with respect to tropho-species was analyzed using network-based null models in econullnetr(16) with the ‘generate_null_net’ command, visually represented with the ‘plot_preferences’ command. Standardized effect sizes of prey choice for each combination of spider genus, sex and life stage, indicative of the extent of deviation from random, were extracted from the null models and compared between genera, sexes and life stages using permutational multivariate analysis of variance (PerMANOVA) via the ‘adonis’ function in vegan(24). To determine any tropho-species-specific differences, these data were further analyzed via similarity percentages analysis (SIMPER), also in vegan.

1.        M. J. Roberts, The Spiders of Great Britain and Ireland (Compact Edition) (Harley Books, Colchester, UK, ed. 3rd, 1993).

2.        H. Krehenwinkel, S. Kennedy, S. Pekár, R. G. Gillespie, 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).

3.        J. P. Cuff, L. E. Drake, M. P. T. G. Tercel, J. E. Stockdale, P. Orozco-terWengel, J. R. Bell, I. P. Vaughan, C. T. Müller, W. O. C. Symondson, Money spider dietary choice in pre- and post-harvest cereal crops using metabarcoding. Ecol. Entomol. 46, 249–261 (2021).

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

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

6.        P. D. Schloss, S. L. Westcott, T. Ryabin, J. R. Hall, M. Hartmann, E. B. Hollister, R. A. Lesniewski, B. B. Oakley, D. H. Parks, C. J. Robinson, J. W. Sahl, B. Stres, G. G. Thallinger, D. J. Van Horn, C. F. Weber, Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).

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

8.        C. Camacho, G. Coulouris, V. Avagyan, N. Ma, J. Papadopoulos, K. Bealer, T. L. Madden, BLAST+: architecture and applications. BMC Bioinformatics. 10, 1–9 (2009).

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

10.      D. H. Huson, S. Beier, I. Flade, A. Górska, M. El-Hadidi, S. Mitra, H. J. Ruscheweyh, R. Tappu, MEGAN Community Edition - interactive exploration and analysis of large-scale microbiome sequencing data. PLoS Comput. Biol. 12, 1–12 (2016).

11.      L. E. Drake, J. P. Cuff, R. E. Young, A. Marchbank, E. A. Chadwick, W. O. C. Symondson, 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.

12.      J. P. Cuff, S. M. Wilder, M. P. T. G. Tercel, R. Hunt, S. Oluwaseun, P. S. Morley, R. A. Badell-Grau, I. P. Vaughan, J. R. Bell, P. Orozco-terWengel, W. O. C. Symondson, C. T. Müller, MEDI: Macronutrient Extraction and Determination from invertebrates, a rapid, cheap and streamlined protocol. Methods Ecol. Evol. 2021, 1–9 (2021).

13.      J. P. Cuff, S. M. Wilder, MEDI: Macronutrient Extraction and Determination from Invertebrates. Protocols.io (2021), p. 49505.

14.      J. P. Cuff, Further micro-scaled MEDI (macronutrient extraction and determination from invertebrates). Protocols.io (2021), , doi:https://dx.doi.org/10.17504/protocols.io.bw5hpg36.

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

16.      I. P. Vaughan, N. J. Gotelli, J. Memmott, C. E. Pearson, G. Woodward, W. O. C. Symondson, 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).

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

18.      N. E. Hamilton, M. Ferry, ggtern: ternary diagrams using ggplot2. J. Stat. Softw. 87, 3 (2018).

19.      H. Wickham, ggplot2: Elegant Graphics for Data Analysis (2016).

20.      G. Brock, V. Pihur, S. Datta, S. Datta, clValid: an R package for cluster validation. J. Stat. Softw. 25, 1–22 (2008).

21.      G. R. Warnes, B. Bolker, L. Bonebakker, R. Gentleman, W. Huber, A. Liaw, T. Lumley, M. Maechler, A. Magnusson, S. Moeller, M. Schwartz, B. Venables, gplots: Various R programming tools for plotting data (2020).

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

23.      S. Garnier, viridis: default color maps from ‘matplotlib’ (2018).

24.      J. Oksanen, F. G. Blanchet, R. Kindt, P. Legendre, P. R. Minchin, R. B. O’Hara, G. L. Simpson, P. Solymos, M. H. H. Stevens, E. Szoecs, H. Wagner, vegan: Community Ecology Package (2016).