Identification of potential insect ecological interactions using a metabarcoding approach

Targeting insect ecological interactions

An ideal method for analysis would apply across a variety of animal targets of interest with minimal methodological variation. While a great many genes have been employed to diagnose species interactions (Pornon et al., 2016; Rowe et al., 2021; Kim et al., 2022), and some primers are designed specifically for this (e.g., Zeale et al., 2011), fewer analyses have compared their use across different systems or hosts. Insects are an ideal model for developing and testing target genes as their great diversity and ecological versatility offer a very wide spectrum of interactions (Yang & Gratton, 2014; Crespo-Pérez et al., 2020). In this work we developed a generalized method for the identification of potential ecological interactions by using multiple target genes for the recovery of interaction data from three insect groups of substantial economic importance in North America: the invasive spongy moth (Lymantria dispar dispar, Lepidoptera) and emerald ash borer (Agrilus planipennis, Coleoptera), as well as a variety of cavity nesting pollinators (Hymenoptera) in their larval form where taxonomic identity is challenging to establish.

Spongy moths

Spongy moths were introduced to North America in 1869 and have caused substantial damage to more than 300 species of trees by defoliation during their caterpillar stage, reducing growth, root, and fruit production (Elkinton & Liebhold, 1990; Liebhold et al., 1995; Davidson, Gottschalk & Johnson, 1999; Muzika & Liebhold, 1999; Kozlowski, Kramer & Pallardy, 1991; Nakajima, 2015). Defoliated trees are often weaker and more susceptible to dying, with both economic and ecological impacts (Davidson, Gottschalk & Johnson, 1999). Defoliation reduces timber production (Leuschner et al., 1996), and in 2021 roughly 1.8 million hectares of forests in Ontario were defoliated by spongy moths (Government of Ontario, 2014). Between timber loss and biological controls, spongy moths cost North America nearly $3.2 billion annually (Bradshaw et al., 2016).

In addition to direct economic costs, spongy moth outbreaks cause disruption of ecological communities and food webs, reducing native moth populations and altering the structure and composition of forests (Gurnell, 1983; Work & McCullough, 2000; Fajvan & Wood, 1996; Canham, 1988). Efforts to control spongy moth populations and to prevent their spread include manually removing and killing the insects, the application of pesticides, and the introduction of biological controls (i.e., natural enemies of the moths, such as the parasitic wasp Ooencyrtus kuvanae and the fungus Entomophaga maimaiga; Mcmanus & Csóka, 2007; Elkinton & Liebhold, 1990). Spongy moth populations fluctuate annually, with large-scale outbreaks occurring synchronously across large distances every 8 to 10 years, and smaller outbreaks every 4 to 5 years (Berryman, 1996; Davidson, Gottschalk & Johnson, 1999; Peltonen et al., 2002; Haynes, Liebhold & Johnson, 2009). The causes of these outbreaks have not yet been determined, but ecological interactions (e.g., with predators and pathogens) have been hypothesized to play a role (Allstadt et al., 2013; Elkinton & Liebhold, 1990). As a result, determining large-scale species interactions may facilitate better predictions of outbreaks and present mechanisms for population control.

Emerald ash borers

The emerald ash borer is a wood-boring beetle native to northeastern Asia. It was first detected in North America near Detroit, Michigan in 2002, although it may have been present undetected since the 1990s (Cappaert et al., 2005; Siegert et al., 2014). Emerald ash borer larvae consume the phloem and cambium of ash trees (Fraxinus spp.) creating S-shaped galleries that disrupt nutrient transport (Cappaert et al., 2005; Wang et al., 2010). Infestations are associated with extensive canopy dieback, especially in green ash (Fraxinus pennsylvanica; Anulewicz, McCullough & Cappaert, 2007), causing tree mortality often exceeding 99% within 6–10 years of beetle arrival (Natural Resources Canada, 2013; Knight, Brown & Long, 2013; Klooster et al., 2014). As a direct result, six North American ash species have been added to the IUCN Red List of Threatened Species (International Union for Conservation of Nature, 2022). Since the beetle’s introduction, tens of millions of trees have been lost (Cappaert et al., 2005) and they are estimated to be the costliest forest pest in North America, causing approximately $850 million in local government expenditures (i.e., tree removal, replacement, and treatment) and $380 million in lost residential property values each year in the USA alone (Aukema et al., 2011).

Ecological impacts include reduced availability of food and shelter for organisms that rely on ash trees, e.g., beavers, rabbits, deer, purple finches, and wood ducks (Schlesinger, 1990), and extinction risks in dependent species, e.g., 43 monophagous native arthropod species like the eastern ash bark beetle (Hylesinus aculeatus) and black-headed ash sawfly (Tethida barda; Gandhi & Herms, 2010a). Like spongy moths, emerald ash borer infestations alter the structure and composition of forests: for example, canopy dieback allows increased light penetration to the subcanopy, encouraging the growth of understory species and altering their distribution and abundance (Canham, 1988; Gandhi & Herms, 2010b; Flower, Knight & Gonzalez-Meler, 2013). Identifying the broad-scale ecological interactions of emerald ash borers may allow for the monitoring of known biological controls or for the identification of novel ones.

Pollinators

Bees are arguably the best-known pollinators, and eusocial bees are heavily managed for efficient, mass agricultural pollination (Kleijn et al., 2015; Klein et al., 2006), due to their domesticity, ease of transport, and large numbers (Iwasaki & Hogendoorn, 2021). As they offer key ecosystem services such as agricultural pollination and, in the case of the honeybee (Apis mellifera), production of honey, they have been a major focus for conservation (Iwasaki & Hogendoorn, 2021). Despite the popularity of these commercial pollinators, most pollinators (bees, wasps, other insects) are not colonial and cannot be managed commercially yet provide crucial pollination services in both natural and agricultural systems (Garibaldi et al., 2011, 2013; Park et al., 2016; Russo et al., 2015). Because they are not colonial, these taxa are more challenging to study. They can be collected in reasonably large numbers as larva but then there is an increased challenge related to their taxonomic identification. Increasing evidence suggests that non-colonial pollinator contributions to agricultural pollination are underestimated and their ecological interactions differ from colonial honeybees requiring alternative management strategies. For example, Gresty et al. (2018) demonstrated that floral provision thought to bolster pollinator populations and mitigate the effects of large-scale agricultural land conversion were not well used by solitary bees, suggesting different ecological interactions. Significant behavioral differences (e.g., shorter flight seasons in native bees; Bosch & Kemp, 2002) and regional differences in community composition (Danforth, Minckley & Neff, 2019) make it difficult to extrapolate ecological data across taxa and region. Analysis of symbiont and holobiont data could provide a rapid mechanism to evaluate ecological interactions where visual observations are not practical, particularly in low density populations with short duration seasonal activity.

In this study we examine potential target regions and protocols to develop a general metabarcoding approach which can be applied to a variety of insects to detect their potential interactions with other taxa. Our objectives are to: (1) evaluate different target regions for the identification of potential ecological interactions in these taxonomic groups, (2) assess whether these analyses should be targeted at different life stages (e.g., eggs, larvae, pupae, and adults), and (3) compare the results of “pooling” polymerase chain reaction (PCR) products from different gene regions into a single sequencing library vs. individual analysis for each marker for analytical efficiency. Direct comparisons of the taxonomic coverage of target gene regions for metabarcoding of potential interactions across a range of insect models will guide methodological choices in ecological analyses.

Sample collection and study design

All samples used in this analysis were acquired from the collections of the Centre for Biodiversity Genomics, University of Guelph, Canada. They had been collected and sorted using clean protocols (gloves, cleaned tools, etc.) to minimize any cross contamination and stored in individual sealed vials. To simplify methodological testing, we developed our protocol using a hierarchical testing design. We started with spongy moths as pooled samples following the approach for Drinkwater et al. (2019). While this does not assess individual interactions, we used this approach to increase DNA yield for initial primer testing. We tested a subset of primers on these samples. We then used an individual approach for a large set of emerald ash borer individuals (where a sample consisted of a whole or partial individual), and larvae from homes deployed to attract native and non-colonial bees, while also reducing reaction volumes and increasing the number of target regions.

DNA extraction

For both spongy moths and emerald ash borers we extracted DNA using the Qiagen blood and tissue kit following the manufacturer’s guidelines, with two modifications. First, specimens were ground using a pestle during initial lysis, and second, we decreased the final DNA elution to 100 µL of Buffer AE for spongy moths and 150 µL for emerald ash borers, to increase DNA yields. For spongy moth extractions we used 4–5 eggs or larvae from the same host tree in order to get sufficient biomass to extract DNA, thus each “sample” is a pool of 4–5 individuals (e.g., see the method used by Drinkwater et al., 2019). For emerald ash borers we extracted DNA from the posterior half of the abdomen for imagoes, the posterior third for pupa specimens, and an approximately 6 mm fragment of the posterior end of larva specimens (or the whole specimen if classified as larva size 1). Extraction blanks were included using sterile water.

Amplification of target gene regions

Sequencing and data analysis

Amplified DNA was sent to the Barts and the London Genome Center at Queen Mary University of London’s Blizzard Institute, with the exception of 18S fragments (see below). All samples sent were processed identically including bead clean up and size selection by 0.9x Ampure Beads (Beckman Coulter, Brea, CA, USA), quantification, quality control (QC) and normalization using Qubit (Invitrogen, Carlsbad, CA, USA) nucleotide quantification, and DNA D100 Tape station (Agilent). All PCRs were independently barcoded using single indexes and sequenced using a MiSeq v3 2 × 300 cycle run (Illumina, San Diego, CA, USA). Raw sequences were demultiplexed on site.

Cutadapt v3.7 (Martin, 2011) was used in paired end mode to separate markers in pooled libraries based on primer pair. Then, these were processed the same way as their individual PCRs correspondents. Demultiplexed ITS and rbcL files were processed using the DADA2 pipeline (Callahan et al., 2016) in R Studio (R Core Team, 2023) following Garrett et al. (2023). We first removed amplicon primers using Cutadapt. In brief, DADA2 was used to merge paired reads, trim sequence length, and remove errors in sequence profiles using the learnErrors function. Chimeras were removed and amplicon sequence variant (ASV) tables were generated for each primer pair and each taxonomic group individually. Resulting ASVs for ITS and rbcL sequences were compared to the nucleotide collection in GenBank using the BLAST algorithm. An identification was retained for further consideration if the identity match was ≥97% (with 100% overlap), and if the detection count was greater than the highest read count in the negative control (negative filtering). We then screened taxonomic identity against available ecological data on distributions (i.e., the organism can be found in the area the samples were collected in according to the Global Biodiversity Information Facility (Global Biodiversity Information Facility (GBIF), 2023). Several potential interactions of specific interests which met some but not all of these conservative criteria are also reported but with cautionary notes (see Results).

16S reads were processed using QIIME2 (Bolyen et al., 2019). Sequences were trimmed using Cutadapt v3.7 (Martin, 2011) and denoised using DADA2 (Callahan et al., 2016). The denoising parameters (trim-left-f, trim-left-r, trunc-len-f, and trunc-len-r) were as follows: 0, 0, 268, and 218 for emerald ash borers; 0, 0, 271, and 202 for spongy moths; 0, 0, 272, and 203 for pollinators. An identification was retained if the detection count was greater than the highest read count in the negative control, if it was not associated with the positive control, and if the taxonomic identity matched known ecological data (i.e., the organism is commonly associated with soil/water/plants or is an insect pathogen according to MicrobeAtlas v1.0 (Matias Rodrigues et al., 2017) and is thus unlikely to be a human/lab contaminant).

Spongy moths and emerald ash borer COI sequences were processed in mBRAVE (Ratnasingham, 2019) using the default parameters except for the following: trim front = 30 bp, trim end = 24 bp, min Quality Value (QV) = 0 qv, max length (pre-trim) = 600 bp, max bases with low QV (<20) = 75%, max bases with ultra-low QV (<10) = 75%, ID distance threshold = 5%, exclude from operational taxonomic unit (OTU) threshold = 5%, read sub-sampling-max reads per sample = 2,500, read sub-sampling-max reads per contig = 200, paired end merging = merge, assembler min overlap = 10 bp, and assembler max substitution = 20 bp. The following system reference libraries were used: Insecta (SYS-CRLINSECTA), Non-Insect Arthropoda (SYS-CRLNONINSECTARTH), Non-Arthropoda Invertebrates (SYS-CRLNONARTHINVERT), and Chordata (SYS-CRLCHORDATA). Pollinator COI was similarly processed but used the CBG Authoritative Canadian Reference Library 2022 (DS-CANREF22) with trim front and end set to 20, ID distance threshold = 1.5%, exclude from OTU threshold = 3%, read sub-sampling − max reads per sample = 2,000, min QV = 10 qv, max length = 1,000 bp, and assembler max substitution = 5 bp. An identification was retained for further evaluation if the mean length fell within the expected range (i.e., expected amplicon length ± 1 bp), if “pc_sim_mean” was ≥97%, and if “id_gaps_mean” equaled zero. We used full negative filtering using the negative controls. All putative taxonomic identifications were compared against known distribution records (e.g., in Global Biodiversity Information Facility (GBIF), 2023). The coleopterans Tylonotus bimaculatus (Cerambicidae) and Hylesinus aculeatus (Curculionidae) were detected in some presumed emerald ash borer samples, and these were removed from further analysis as larvae were likely misidentified during the collection stage. For each larval sample, the ASV with the highest read count above 100 reads was considered the “occupant” of the sample (Handler et al., 2024). If all read counts were lower than 100, or there were two very similar IDs it was labelled as an “unclear ID” sample reflecting poor amplification or ID in general. All other associated ASVs were considered potential interactions. The most common occupants (found in a minimum of three samples) were chosen as a focus for analysis.

Because the 18S amplicons exceed the length which can be processed on Illumina platforms (>900 bp for the target nematodes), 18S emerald ash borer amplicons were sequenced using the MinION platform (Oxford Nanopore, Oxford, UK) with library prep carried out as described in Floyd, Prosser & Jafarpour (2023) and run on a Flongle flow cell. Data was analyzed using the methods detailed in Hebert et al. (2025). To assign putative taxonomic identification to sequences, the SINTAX tool (Edgar, 2016) was employed, with a custom reference library containing the nematode sequences from 18S NemaBase (Gattoni et al., 2023) and sequences of all other taxa from the SILVA public database.

All figures were generated in R (R Core Team, 2023) using the ggplot2 package (Wickham, 2016) and arranged into multi-panel figures in power point.

Spongy moths

We processed 1,460,234 COI, 199,979 ITS, 539,627 rbcL and 292,782 16S raw reads. These were reduced to 281 non-host reads (spongy moth data was excluded), 92,616, 377,885 and 18,091 filtered reads which were converted to 15, 167, 42 and 97 BINS (COI) or ASVs respectively. We identified 126 taxa across the four markers (ITS = 49; 16S = 41; COI = 8; rbcL = 28). Of these, 119 were unique taxa (plants identified by both rbcL and ITS were not counted twice). Most identifications were made to level of genus, and these represented 60 different orders, with two unclassified algae and an uncultured microbe included in this count. We identified 19 orders (six fungi and 13 plants) using the ITS marker, detecting 17 in the eggs and 14 in the larvae (Table S1). Using the 16S marker we identified 28 bacteria orders, with 22 and 20 orders detected in the egg and larvae life stages, respectively (Table S2). We detected Wolbachia sp. (order Rickettsiales), an insect symbiont, in five samples. Using the COI marker, we detected six orders of insects, other than the host itself, with five orders in the egg stage and four in the larval stage (Table S3), which we suspect represent potential parasites. We detected 15 plant orders with rbcL (eight of which were also detected using ITS), with eggs containing plant DNA from 13 orders, and larval samples containing plant DNA from 11 (Table S4).

These detections from 60 orders accounted for 317 actual occurrences across the two life stages (Fig. 1). There were 167 order-level occurrences in the egg life stage (nine animal, 60 bacteria, 21 fungi, and 77 plant), while we recorded 150 occurrences in the egg life stage (six animal, 64 bacteria, 14 fungi, and 66 plant).

Emerald ash borers

We processed 11,361,252 COI, 2,443,232 ITS, 6,173,341 16S, and 287,180 18S raw reads. These were reduced to 114,088 non-host reads (emerald ash borer data was excluded), 942,740, 8,154 and 8,126 filtered reads which were converted to 6, 36, 38, and 31 BINS (COI, 18S) or ASVs respectively. We identified 65 taxa across the four markers (ITS = 23; 16S = 21; COI = 4; 18S = 17). These represented 61 unique identifications from 30 orders. We detected 11 orders of fungi using the ITS marker, with larva size 1 showing a richness of five orders, larva size 2 seven orders, larva size 3 three orders, pupae six orders, and imagoes eight orders (Table S5). Notably, we detected the fungal biocontrol agent Beauveria sp. (order Hypocreales) in two samples. However, it should be noted that this was also detected in one of the negative controls. We retained the occurrence here because it was present in very small amounts in this control compared those detected in the samples and there is no other known source of the fungus in our processing. It was also independently verified by 18S data (see below) supporting it as a true positive. Using the 16S marker, we detected 13 orders of bacteria. We identified five different orders in the larva size 1 stage, three in the larva size 2 stage, two in the larva size 3 stage, eight in the pupal stage, and nine in the adult stage (Table S6). We detected the insect symbiont Wolbachia sp. (order Rickettsiales) in one sample. We also detected four arthropod orders using the COI marker, including potential parasites. Two orders were detected in the larva size 2, one in the larva size 3 and pupal stage, and two in the adult stage. Notable identifications included Ichneumonidae parasitoid wasps (Table S7). With the 18S marker, nine orders were detected, with six of the nine representing fungi (all overlapping with those identified using ITS), and the remaining three representing insects, nematodes, and a parasitic alveolate protist (Cryptosporidium sp., order Cryptosporida): three in the larva size 1 stage, five in the larva size 2, one in the larva size 3, four in the pupae, and two in the imagoes (Table S8).

These detections from 30 orders accounted for 259 actual detections across life stages (Fig. 2). There were 29 potential order level interactions detected in larva size 1 (three animal, 17 bacteria, and nine fungi), 64 in larva size 2 (eight animal, 24 bacteria, and 32 fungi), 10 in larva size 3 (two animal, three bacteria, and five fungi), 74 in the pupae (five animal, 31 bacteria, and 38 fungi), and 82 in the adult beetles (four animal, 48 bacteria, 29 fungi, and one parasitic alveolate).

Pollinators

We processed 2,054,675 COI, 1,395,445 ITS, 1,649,988 rbcL, and 2,731,151 16S raw reads. These were reduced to 221,846, 1,088,989, 1,564,079 and 2,220,678 filtered reads which were converted to 188, 1,538, 727 and 3,169 BINS (COI) or ASVs respectively. We identified 330 taxa from 89 orders across the four markers (rbcL = 22; ITS = 24 plants and 15 fungi; 16S = 20; COI = 8 Tables S9–S12). These represented 73 unique order-level identifications. A total of 16 orders of plants were identified by both rbcL and ITS and were not counted twice. We identified 24 species of Hymenoptera among the larvae, the most common of which were the bees Heriades carinatus, Megachile campanulae, M. centuncularis, M. relativa, Osmia tersula, and the wasp Symmorphus bifasciatus. Other COI identifications included known parasites and prey items (e.g., of the wasp Apocrita sp. which also uses solitary bee homes). We identified 17 fungal genera of which the most common was Alternaria, followed by Ascosphaera, Aspergillus, and Penicillium. There was overlap between plants identified using ITS and rbcL but much higher richness estimates from ITS (e.g., all individuals had more than 20 genera of plant identified with ITS while this was a maximum identified with rbcL). The most common plant genera identified were associated with the Asteraceae, Brassicaceae, Fabaceae, and Rosaceae. Using the 16S region we identified 39 bacterial genera, the most common being Wolbachia, followed by Apilactobacillus and Mellisococcus.

Individual vs. pooled samples from pollinators

Given that the same PCRs were used for pooled and individual libraries, the similarity was surprisingly small between them. Taxonomic coverage for pooled PCRs was much lower than individual PCRs in all cases (Fig. 3). The best coverage for pooled PCRs was offered by the COI marker, where 15 genera were detected in pooled PCRs vs. 24 in the same samples when not pooled. Comparatively, the other markers had more notable differences in coverage. Individual ITS plant and fungal amplifications identified 180 genera but only 45 were detected in the pooled samples with one taxa, Parthenocissus unique to pooled data. Within the 16S data, two genera were identified in pooled samples compared to the 39 identified in individual samples. Pooled COI data recovered far fewer host larval IDs, with most not meeting the criteria for identification of a host. The identified genera for the pooled samples frequently did not include most common genera in the individual samples.

Confirming target insect identification

DNA barcoding is an efficient and well-known method for the identification of insects (Jung, Duwal & Lee, 2011; Park et al., 2011; Boehme, Amendt & Zehner, 2012), and COI is widely used as the primary marker for taxonomic identification of animals (Hebert et al., 2003). Here we use it with a dual purpose: to confirm the identity of the host and to identify potential metazoan parasites. Unsurprisingly, host DNA dominated the sequencing datasets, with only trace amounts of DNA from other organisms, in terms of relative read abundances (RRA). We did use strict filtering protocols from negative controls which might have resulted in artificially reduced detections of non-host taxa. The alternative, a more relaxed filtering, or not filtering but tracking potential contaminants, could raise issues around non-relevant detections (see below discussion of environmental DNA). Among emerald ash borer samples, we detected four orders of arthropods other than the host itself: two families of mites (Oppiidae, order Sarcoptiformes, and Eupodidae, order Trombidiformes), and one family of gall midges (Cecidomyiidae, order Diptera). Gall midges are primarily plant pests, and while some species can act as parasites of other arthropod larvae (Hayon, Mendel & Dorchin, 2016), we could not identify these further. These are likely co-occurring taxa rather than parasites. Notably, we detected a family of parasitoid wasps (Ichneumonidae, order Hymenoptera) in larva size 2, larva size 3, and pupa. Ichneumonids have previously been observed to parasitize emerald ash borers (Duan et al., 2009) making these highly likely cases of true parasitism.

The widest variety of taxa was found in the pollinator data where 24 occupants were detected, including aphids (probably kleptoparasites), and the wasp Apocrita sp., which may use the homes set out for non-colonial bees. The criteria used to determine what constitutes an “occupant” was arbitrary, and read abundances may vary depending on many factors other than true abundance and biomass (such as sequencing depth, primer-binding efficiency, etc.). However the genomic extracts used here had previously been screened and were selected based on having a clear “signal”, thus we hoped to limit the challenge. However, we still encountered a small number where the ID of the occupant was unclear. In part this is due to a much larger reference library being available in this analysis than in a previous one. It may also be due to a common challenge where hymenopteran DNA is notoriously hard to sequence when in mixed templates. This is hypothesized to reflect the extremely high AT content of the mitochondrial genomes of Hymenoptera (Clare et al., 2008) which may lead to weaker PCR amplification in mixed templates where a more GC rich template is available. Among those where ID was clear, we did recover a reasonable set of taxonomic identifications, with common home occupants correctly labeled as “occupant”, and known parasites as secondary signals. We suggest that these criteria be adapted to case-specific values for each new analysis as any such criteria will vary with sequence depth, primer binding etc., and analysis for ecological objectives would be better aimed at screening unpooled larva, even though this increases analytical costs, to better control for occupancy by multiple species.

The detection of potential insect-plant interactions

We used two markers to detect plant taxa: ITS and rbcL. ITS was amplified using two different primer pairs in our study. While the forward primer ITS3 was designed specifically for fungal taxa (White et al., 1990), the primer pair ITS-S2F/ITS4 is known to cross amplify fungi and plants, making it a useful general marker (Chen et al., 2010; Cheng et al., 2016). It may be most appropriate for use when plants are not the only target. For example, we used this primer set with spongy moths where there were a variety of host trees, and we were also interested in fungal amplification. In this case we did amplify a variety of plants and fungi, but the host tree was only identified once. While we were not specifically targeting plants, it was somewhat surprising that the host was not often recovered, and may be linked to the strict filtering of our data, and the age and preservation of the specimens (extracted more than two years after collection and specimens not frozen immediately). In contrast to ITS, rbcL does not readily cross amplify most taxa outside of the plants, other than green algae or cyanobacteria (neither of which were expected in our samples). In spongy moths, rbcL identified 15 orders of plants compared to 13 with ITS, but only eight orders were common between the two. Unlike ITS, rbcL did successfully detect the host tree in 11 of the 18 samples. In pollinators the number of detections was much higher with ITS, recovering far more taxonomic diversity than rbcL. Individual pollinator larval samples contained DNA identified as up to 20 plant genera (rbcL) or more than 20 plant genera (ITS).

The detection of potential insect-fungal interactions

The primer pair ITS3/ITS4 is less likely to amplify plant taxa and was used to amplify fungal DNA from emerald ash borers, where the tree host (ash) was known and, thus, devoting sequencing power to plant identification was not useful. Additionally, ITS3/ITS4 primers are known to have better specificity and amplify a greater variety of fungal taxa in comparison to other ITS primers (Yu et al., 2022). Most of the fungi detected within the emerald ash borer samples are associated with soil (e.g., Fusarium sp.), decaying plant matter (e.g., Yamadazyma sp.), lichens (e.g., Candelaria sp.), or plant pathogens (e.g., Diaporthe sp.), and not likely the larvae, pupae, or imagoes themselves, suggesting this is co-occurrence data rather than direct interactions, though it may be ingested material. One of the most commonly identified fungi is Alternaria sp., a genus of plant pathogens, involved in decomposition, but also ubiquitous in the environment and found almost everywhere (and considered a general environmental contaminant of airborne eDNA studies, E. Clare personal observation). However, one fungus of interest that was detected in two samples (one larva size 2 and one adult beetle) is Beauveria sp. (order Hypocreales). Beauveria spp. are entomopathogenic fungi that are known to successfully infect and neutralize emerald ash borers. Although Beauveria occurs naturally in the environment (e.g., in soil; Rehner et al., 2011), biocontrol products containing Beauveria (e.g., FraxiProtec) are also commercially available and have been shown to significantly reduce emerald ash borer populations (Srei et al., 2020). This detection here suggests some biocontrol may be ongoing and detectable by this analysis.

The detection of potential insect-microbial interactions

Similar to the fungi detected by ITS, most of the bacteria associated with our samples are ubiquitous throughout the environment, and can be associated with soil (e.g., Massilia sp. in spongy moth samples, and Terriglobus sp. in emerald ash borer samples), plants (e.g., Capsulimonas sp.), lichen (e.g., Lichenihabitans sp.), or the rhizosphere (e.g., Neorhizobium sp. found in emerald ash borer samples). However, one notable bacterium found in all hosts is Wolbachia sp. (order Rickettsiales). Wolbachia sp. are intracellular parasites that can infect many insects and nematodes (Tagami & Miura, 2004), and are associated with reproductive abnormalities, such as the feminization of genetic males, parthenogenesis induction, and reproductive incompatibility (i.e., through cytoplasmic incompatibility between hosts infected by different strains or between uninfected/infected hosts; reviewed in Werren, 1997). However, it has recently been shown that Wolbachia might confer some protection to their host against viral infections, at least in dipterans (Cogni et al., 2021; Pimentel et al., 2021). Wolbachia’s ability to induce genetic incompatibility in populations coupled with the ability for hosts to transmit the infection vertically and horizontally means it may have potential biocontrol applications (Werren, 1997; Zabalou et al., 2004).

Challenges of detecting parasites from hosts

We attempted to use both COI and 18S to detect parasitism. Parasite detection is extremely challenging for a number of reasons. First, differentiating true parasites from environmental DNA contamination of the body is nearly impossible. Second, without destruction of the sample, detecting internal parasites is rare (e.g., sometimes detected from faeces or during voucher saving extraction methods) and most signals will likely be from surface contamination. A challenge with COI is the swamping of secondary signals with host DNA which is both more abundant and co-amplified, and most primers in DNA barcoding have been designed to minimize taxonomic bias. An alternative target is 18S where some primers show higher affinity for nematodes, intracellular parasites, etc. However, short 18S regions suitable for most high throughput sequencing on platforms such as Illumina provide limited taxonomic resolution. To address this gap, we sequenced long 18S amplicons on the MinION instrument (Oxford Nanopore, Oxford, UK) which provides much lower sequencing depth but can accommodate long read lengths at low cost, even with a high negative detection rate.

We tested this approach on emerald ash borer specimens and, interestingly, the 18S region detected the greatest diversity of interaction types (Fig. 2), including animals (nematodes and arthropods, which were expected), but also some fungi, and a parasitic alveolate. All fungal orders detected using 18S were also detected using ITS, but, despite being a longer amplicon length, the taxonomic resolution offered by 18S was the same or worse than ITS, likely due to a combination of generalist primers lacking resolution and incomplete reference libraries (e.g., taxa recovered by both markers but assigned with ITS to genera such as Diaporthe sp., Niesslia sp., and Cryptococcus sp., were only resolved at a higher taxonomic level with 18S). This is consistent with other studies that have found the ITS marker to provide a more extensive view of fungal diversity and allow more precise identifications to be made compared to the 18S marker (Lord et al., 2002; Liu et al., 2015). One notable detection was the family Cordycipitaceae (order Hypocreales), which contains the previously noted fungal biocontrol genus Beauveria, independently confirming its presence in the same samples, and providing support to its detection when we treated it with caution in ITS data. While the recovery was much lower in taxonomic richness, the lower resolution of the 18S region might be a useful screening tool for other interactions.

The non-host arthropods detected belong to order Hymenoptera. Potential interactions with nematodes (order Rhabditida) detected include Rhabditolaimus sp. (bacteriophagous and commensal with other beetles, e.g., Scolytus multistriatus; Ryss & Polyanina, 2022), Panagrellus sp. (free-living nematodes associated with many habitats, including insect frass; Ferris, 2009; Srinivasan et al., 2013) and Neotylenchidae (members can be parasitic, e.g., Deladenus proximus; Zieman et al., 2015). Finally, the parasitic alveolate detected (Cryptosporidium sp.) has been reported to infect humans and animals including insects (Helmy & Hafez, 2022).

Novel challenges in detecting species interactions

Environmental DNA (eDNA) refers to DNA shed by organisms (e.g., through cells, excreta, etc.) that is deposited throughout the environment, including in soil (Marquina et al., 2019; Foucher et al., 2020), water (Uchida et al., 2020), and air (Clare et al., 2021; Clare et al., 2022). eDNA can also be spread between organisms directly through contact. For example, insects leave eDNA behind on the surface of plants they inhabit/consume (Kudoh, Minamoto & Yamamoto, 2020; Allen et al., 2023) and this can be picked up by another insect without any actual interaction. This was demonstrated by Huszarik et al. (2023), who found insect eDNA was spread to spiders who occupied the same wet pitfall trap as an exotic insect that was not consumed by the spiders. It is thus inevitable for target organisms to become coated in eDNA simply by occupying their natural environment or during collection/storage depending on the medium (e.g., pitfall traps, malaise traps, tubes, ethanol, etc.; Shokralla, Singer & Hajibabaei, 2010; Huszarik et al., 2023). While eDNA’s ubiquitous nature allows for its use in biodiversity and community composition assessments (Clare et al., 2022; Macher et al., 2023), it can complicate analyses of interactions through metabarcoding. More specifically, it is not clear how long eDNA remains detectable and it is unclear whether eDNA contamination on the surface of insects or other target organisms poses a challenge to the assessment of species interactions. Surface contamination could make it difficult to determine true interactions from mere co-occurrence, but only if eDNA regularly persists, something which has not been robustly tested.

Some of the detections recorded here were unexpected. For example, when attempting to confirm target hosts, we often detected other insects which are not likely directly associated with the target. In spongy moth samples we detected six insect orders, including other plant pests (e.g., Epuraea sp. and Pineus sp.) and predatory insects (e.g., Hemerobius sp. and Lestodiplosis sp.). Spongy moths may encounter other plant pests (either directly or indirectly) on the host trees and may have even been predated on by some of the insects detected like lacewings (Hemerobius sp., order Neuroptera). For example, green lacewings are generalist predators that have been observed to feed on lepidopteran eggs and small larvae (Brown & Cameron, 1982; Tauber, Tauber & Albuquerque, 2009). However, it is impossible to confirm whether these are predatory interactions or simply encounters with other insects in the environment. Interestingly, in one spongy moth sample from a pine tree a pine-feeding aphid (Pineus sp.) and its predator (Lestodiplosis sp.) were also detected (Gagn & Havill, 2020). While it is unlikely this was a direct interaction with the spongy moth, it suggests the detection of another local interaction occurring on the host tree, a fascinating ecological process to be recorded in surface contamination by environmental DNA of nearby organisms, though unrelated to the objective of the analysis. Despite the potential accumulation of DNA signatures with time from environmental sources, the actual profile of co-occurring DNA signatures was remarkably stable with life stage (e.g., see Fig. 4 for spongy moths) which suggests this surface contamination risk may be minimal.

The effect of potential eDNA contamination was larger in plant identifications, and this may reflect the more robust cells and pollen of plants. We commonly identified DNA from trees and plants that were not noted during spongy moth collection with both ITS (e.g., Salix sp., Poa sp., etc.) and rbcL (e.g., Rosa sp., Picea sp., etc.). These are possibly a result of pollen from nearby plants sticking to the fuzzy surface of egg masses or to the hairs/surface of the larvae and, in the case of the larvae specifically, we may also be detecting trees that previously hosted the larvae before collection. Among bees it may be impossible to differentiate direct provision of the larva by the mother, material used in nest construction, and casual environmental contact. When the concept is expanded to colonial honeybees, the honey itself is a well-known reservoir of plant and microbial DNA (Hawkins et al., 2015; Balzan et al., 2020; Ribani et al., 2020; Pathiraja et al., 2023).

Most of the fungi detected are ubiquitous throughout the environment and are often found in soil (e.g., Sporobolomyces sp.) or plants (e.g., plant pathogens like Alternaria sp.). However, other fungi (e.g., Punctelia sp.) and algae (e.g., Coccomyxa sp., Trebouxia sp., and Symbiochloris sp.) that often form lichen complexes (Sanders & Masumoto, 2021) were also identified as part of the spongy moth samples by ITS, rbcL, or both. Thus, many identifications made from whole bodies are not target interactions but seem to reflect surface level contamination from co-occurring species through environmental DNA.

There are several potential solutions to this problem of surface contamination. First, it may be prudent to wash specimens prior to DNA extraction (i.e., in a bleach solution) to remove surface DNA (Binetruy et al., 2019; Hausmann et al., 2021; Huszarik et al., 2023), particularly if they have been stored in a common location (e.g., a malaise trap collection pot). This was done to the Hymenoptera larva prior to their original extraction and is thought to remove surface DNA but may limit the detection of fungi and microbes associated with that ecology, though we noted non-target identifications there as well. Similarly, the gut of an insect might be dissected out to determine ingested plants and gut microbiota and parasitism, which are more likely to reflect direct interactions than the causal contact detected by surface environmental DNA. This would also likely decrease the false negative rate in the detection of gut parasites like nematodes by better exposing them to DNA extraction and reducing the proportion of host DNA in the mix. This, however, would cause substantial or complete destruction of the specimen. In marine systems, where diet analysis must consider actual ingested material vs. DNA in swallowed water, one suggested solution has been to analyze DNA in water samples along with DNA from stomach contents to evaluate the potential errors (William. O.C. Symondson, 2012, personal communications). Some version of this, considering body surface and gut separately, bee home tube and larva, or host tree tissue separately from individual pest might similarly help establish likely error rates, though this would also increase the false negative rate: if individual A interacts with parasite B, both may leave DNA in the environment. Thus, removing signals that are present in the environment may be overly conservative. Similarly, restricting our detections only to those for which an ecology has already been established limits the discovery of new true interactions. It may instead be argued that the definitions of symbiome vs. holobiont, interactions vs. co-occurrence are too strict a categorization of relationships which likely fall along a continuum of effects. Strongly interacting species, those with clear coevolutionary “Red Queen” relationships fall at one end, and neutral co-occurring species at the other, but between them are a range of weak interactions, secondary interactions, and tertiary interactions, which form the more complex roles within ecosystems.

Pooling of individual PCRs

One consideration in multi-gene analyses is the cost associated with tracking individual samples with unique tags during sequencing (frequently called indexing). Pooling target regions before indexing samples significantly reduces the costs of analysis, but there are challenges in PCR 2 with the addition of such indexes (e.g., shorter amplicons can become overrepresented during PCR 2). We tested this hypothesis by comparing pooled and unpooled samples from our most taxonomically complex PCRs generated from larval pollinators. Our data suggest that pooling greatly reduced taxonomic representation with many of the most common taxa in unpooled samples completely missing in pooled ones. We observed the same effect in all gene regions suggesting it was not confined to underrepresentation of specific amplicons. The most likely explanation is that pooled libraries lacked sequencing depth compared to individual-PCR libraries: while the pooled libraries contained four times more information (the four markers used for each sample), they were added equimolarly to the sequencing mix, meaning that each pooled library was allocated only one fourth of the sequencing power (which could translate to lower “molecular sampling effort”) compared to individual-PCR counterparts. A more efficient approach may be to amplify samples using PCR 1 primers with unique codes or tags added at the 5′ or 3′ end so that they can be indexed with the same tag during PCR 2 and pooled after (Binladen et al., 2007). This creates a bioinformatic challenge but would reduce indexing steps, albeit increasing laboratory costs (many versions of the same primer, one with each tag). Alternatively, one-step PCR approaches where fusion primers (which include the primer sequence, a sequencing adapter, and nucleotide tags) are used can also eliminate the need for indexing because of the presence of nucleotide tags at the cost of reduced PCR efficiency and increased primer cost (Bohmann et al., 2022). Multiplexing approaches where multiple primers are used to amplify different target regions simultaneously can also be employed to increase efficiency (Zhang et al., 2018). However, difficulties may arise in multiplexed PCRs, such as the formation of primer-dimers and the over-representation of some amplicons over others (Khodakov, Wang & Zhang, 2016), particularly where different extension times are required.