Unveiling invasive insect threats to plant biodiversity: Leveraging eDNA metabarcoding and saturated salt trap solutions for biosurveillance

RESEARCH ARTICLE Unveiling invasive insect threats to plant biodiversity: Leveraging eDNA metabarcoding and saturated salt trap solutions for biosurveillance Yoamel Milia´ n-Garcı´aID1,2*, Cassandre Pyne1, Kate Lindsay1, Adriana Romero1, Robert H. Hanner1 1 Department of Integrative Biology, University of Guelph, Guelph, Ontario, Canada, 2 Ecological and Regulatory (Ecoreg) Solutions Inc., Guelph, Guelph, Ontario, Canada * ymiliang@uoguelph.ca, yoamelmg@gmail.com Abstract The negative global impacts of invasive alien species (IAS) on biodiversity are second only to habitat loss. eDNA metabarcoding allows for a faster and more comprehensive evaluation of community species composition, with a higher taxonomic resolution and less taxonomic expertise required than traditional morphological-based biosurveillance. These advantages have positioned eDNA metabarcoding as the standard method for molecular-based detec- tion of invasive alien species, where fast and accurate detectability allows prompt responses to mitigate their adverse effects. Here, eDNA metabarcoding is used for biosurveillance of invasive alien species regulated by Canada in high-risk areas with four main objectives: i) validate the effectiveness of eDNA metabarcoding of salt trap solutions as a molecular tech- nique for IAS detection, ii) compare detection from DNA extracts obtained from filter quar- ters versus whole filters, iii) benchmark two different bioinformatic pipelines (MetaWorks and mBRAVE), and iv) compare canopy and ground level trapping. eDNA from up to five

IAS (Agrilus planipennis, Daktulosphaira vitifoliae, Lymantria dispar, Popillia japonica, and Trichoferus campestris) were successfully detected across years from 2017 to 2022 in southern Ontario, Canada, with successful morphological validation for all except Lymantria dispar and Trichoferus campestris. Analysis of filter quarters in contrast to whole filters was demonstrated to be insufficient for effective IAS detection in each sample. All IAS were detected in only one filter quarter, suggesting a patchy eDNA distribution on the filter. The MetaWorks and mBRAVE bioinformatics pipelines proved effective in identifying IAS, with MetaWorks yielding a higher success rate when comparing molecular and morphological identifications. Ground-level and canopy-level sampling showed differential IAS recovery rates based on the molecular detection, which also varied per collection year, with all found IAS detected at the canopy level in 2022 while only one (Lymantria dispar) in 2020. The present study ratifies the efficacy and importance of eDNA-based detection in a regulatory context and the utility of adding eDNA metabarcoding of saturated salt trap solutions, a criti- cal tool for IAS detection. PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 1 / 24 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Milia´n-Garcı´a Y, Pyne C, Lindsay K, Romero A, Hanner RH (2023) Unveiling invasive insect threats to plant biodiversity: Leveraging eDNA metabarcoding and saturated salt trap

solutions for biosurveillance. PLoS ONE 18(8): e0290036. https://doi.org/10.1371/journal. pone.0290036 Editor: Ruslan Kalendar, University of Helsinki: Helsingin Yliopisto, FINLAND Received: May 8, 2023 Accepted: July 31, 2023 Published: August 11, 2023 Copyright: © 2023 Milia´n-Garcı´a et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: The raw data (FASTQ files) generated or analyzed during the current study was made publicly available in the Borealis Research Data Repository (2017 collection: https:// doi.org/10.5683/SP2/UM6RYF, 2018 collection: https://doi.org/10.5683/SP2/HTSZNS, and 2019, 2020, and 2022 collections: https://doi.org/10. 5683/SP3/AGUXKR). All images of the IAS specimens discovered in the 2022 collection have been deposited in the Barcode of Life Data Systems (BOLD): dx.doi.org/10.5883/DS-IAS2022.

Introduction The negative effects of invasive alien species (IAS) are among the top human-mediated drivers of biodiversity loss [1]. Environmental DNA (eDNA) metabarcoding has become the method of choice for molecular-based biosurveillance [2, 3]. It allows for cost-effective, multi-species identification in diverse communities by sampling the DNA shed by the organisms into their environment [4]. Compared to conventional morphology-based identification methods, eDNA metabarcoding can reduce the time spent on species identification, sample sorting and processing, and increase the target taxa’s scope by utilizing High-Throughput Sequencing (HTS) technologies coupled with DNA reference databases. The existence of DNA reference databases (e.g., The Barcode of Life Data Systems [BOLD] [5]), which links genetic and mor- phological information, reduces the necessity of taxonomic expertise by assigning known iden- tities to unclassified organisms based on molecular data. eDNA metabarcoding has the capability to detect rare or elusive species from complex envi- ronmental samples, including species that may be considered pests, vectors of diseases, or spe- cies of conservation concern [6]. Additionally, eDNA metabarcoding causes less ecosystem disturbance than traditional biosurveillance methods and can potentially maintain the mor- phological integrity of detected organisms [7, 8]. The non-destructive capacity of eDNA sam- pling is critical in a regulatory context; it allows for specimens to be preserved for secondary

morphological confirmation, which is crucial for the validation of eDNA-methods through comparison to morphology-based identification biosurveillance. Additionally, taxonomic assessment of immature stages of most arthropods is challenging; characters for identification are often dependent on adult morphology and reproductive structures [9]. DNA-based molec- ular identification tools such as eDNA metabarcoding provide a valuable resource for rapidly and accurately identifying morphologically indistinct IAS [10], as eDNA-based classification is independent of morphological integrity and remains invariable regardless of organism’s life stage. Biological invasions are responsible for tremendous global impacts, including substantial biodiversity and economic losses, and management expenditures [11]. Invasive insects cause long-term environmental effects, including reduced ecosystem stability and loss of native spe- cies [12]. Invasive alien species (IAS) can destroy about 400,000 ha of forest annually in Can- ada [13] and are considered second only to wildfires regarding adverse effects on forests. Previous research has shown that invasions of insects alone cost a minimum of US$70.0 billion per year globally [11]. In Canada, the government is committed to strengthening the regulatory framework to control and prevent IAS. The Canadian Food Inspection Agency (CFIA) is primarily responsi- ble for the biosurveillance of terrestrial IAS [14]. Canada has at least 272 federally terrestrial regulated species, including all the potentially invasive species, where the vast majority are

insects [N = 83] (https://inspection.canada.ca/plant-health/invasive-species/regulated-pests/ eng/1363317115207/1363317187811). Two (Agrilus planipennis and Lymantria dispar) out of those 83 insects are among the top IAS that cost the province of Ontario the most money [15]. Within this framework, IAS DNA-based detection, and the need to integrate molecular-based technology into regulatory policy, are continuously being explored to establish cost-effective protocols for the molecular biosurveillance [4, 16]. The main limitation for faster and more routine integration of DNA-based surveillance into current IAS biosurveillance toolkits has been linked to the lack of a thorough technological validation [17]. Four main validation cate- gories have been recognized to support the efficacy of eDNA-based methods: experimental val- idation, comparative morphological/ecological validation, bioinformatic validation, and statistical validation [17]. In the present case study, we aim to address three of these categories. PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 2 / 24 Funding: UoG-BIO acknowledges the financial support of the Canadian Food Inspection Agency (CFIA). [FAP#2122-002 ] The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

The first is to validate eDNA findings by comparing results from filter quarters versus whole filters and exploring whether the former resembles the same detection rate as the latter (experi- mental validation). Second, to compare molecular-based detection to traditional morphologi- cal-based classification (comparative morphological/ecological validation). The third is to compare outputs generated from two informatics pipelines for DNA data processing (bioin- formatic validation). Validation of bioinformatics pipelines can be daunting when considering the plethora of options available for DNA metabarcoding data analyses [18]. Among other things, bioinformatic pipeline selection is linked to the users’ needs and informatics expertise [18]. In a CFIA regulatory context, user-friendly pipelines can support end-user eDNA valida- tion by officially recognizing its utility and acceptable standards for case-specific studies such as IAS detection. In this context, we chose two bioinformatic pipelines (MetaWorks [19] and mBRAVE [https://www.mbrave.net/] [20]) out of the more than thirty-one available [18]. The Multiplex Barcode Research and Visualization Environment (mBRAVE) is a multi-user and friendly platform that supports storing, validating, and analyzing metabarcoding data. This free server-based system is built on the BOLD platform and supports rapid species identifica- tion without demanding a specific user’s informatics infrastructure or expertise. It also offers a secure collaboration and sharing environment [21], which can be essential for end-user eDNA

validity recognition, especially for industry collaborators and regulatory agencies. On the other hand, MetaWorks is a free and open-source software for scalable metabarcoding data analyses of the most popular DNA barcodes [19]. It uses a naive Bayesian classifier and allows sequence composition-based taxonomic assignments [19] with associated confidence mea- sures (bootstrap support values) [19, 22] but also requires a basic bioinformatics background and infrastructure. In the present study, eDNA metabarcoding was used for the molecular biosurveillance of Canadian-regulated terrestrial IAS in high-risk areas with four main objectives: i) validate the effectiveness of eDNA metabarcoding as a molecular technique for IAS detection during a multi-year sampling effort conducted by the CFIA within a regulatory context, ii) compare the capacity of IAS detection when using DNA extracts obtained from filter quarters versus whole filters; iii) benchmark two different bioinformatic pipelines for metabarcoding data analysis (MetaWorks [19] and mBRAVE [https://www.mbrave.net/] [20]) within the context of IAS detection; and iv) compare the capacity of IAS detection when sampling at the canopy and ground levels. To reach these goals, the present study employs an eDNA metabarcoding approach on saturated salt trap solutions (collection fluid) collected from Lindgren funnel traps and focuses on detecting IAS regulated by Canada, adopting a method similar to the one used in previous studies [4, 23]. Materials and methods Target species

The CFIA monitors species considered both invasive and regulated pests. The target species for the current study included all invasive insects regulated by the CFIA (Accessed on May 25, 2022). In addition to invasive species, the list of pests regulated by the CFIA includes 29 fami- lies and 64 genera of insects, five families and seven genera of bacteria, 22 families and 31 gen- era of fungi and 13 families and 23 genera of plants. For simplicity, we here refer to both CFIA-regulated invasive species and regulated pests as “invasive alien species” (IAS). The pres- ent study targeted wood-boring beetles considered as IAS, including long-horned beetles (Cer- ambycidae) (e.g., Anoplophora glabripennis, Anoplophora chinensis, Tetropium fuscum, Tetropium castaneum, and Aromia bungii). Additional species were added to the regulated species list based on noted exceptions in the CFIA list. For example, the CFIA-regulated PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 3 / 24

species list states that Phytophthora species are regulated except for a listed subset. Therefore, the additional regulated but not listed species were added (see S1 Table). Collection locations Lindgren funnel traps were deployed at four undisclosed locations in southern Ontario, Can- ada, in line with the CFIA’s 2019 regulatory survey program aimed at detecting insects intro- duced to Canada through high-risk pathways. The selected locations were Oshawa, Chatham- Kent, Brantford, and Brampton. For each location, six traps were strategically placed at six dif- ferent sites near tree species known to be hosts of the target invasive alien species (IAS). This process was repeated twice, resulting in a total of 48 field samples collected in 2019. In the sub- sequent year, Lindgren funnel traps were set up in two additional southern Ontario locations: Greater Toronto Area and Hamilton. Similarly, six traps were placed at two sites on two differ- ent dates, yielding 24 samples. Furthermore, three more locations (Peterborough, Windsor, and Fort Erie) were sampled twice in 2022, resulting in 30 additional samples. To enhance the analysis, datasets from traps collected in 2017 and 2018 were also included [4, 23]. S1 Fig. con- tains an interactive map displaying the distribution of all the collected sites from 2017 to 2022. Please refer to the supporting material for further details. Collection protocol Lindgren funnel traps were set up with Ultra High Release (UHR) Ethanol Lures from Synergy Semiochemical Corp, including fuscomol and fuscomol acetate as baits. To preserve trapped

insect specimens and samples for eDNA analysis, the collection cups were filled with a satu- rated saltwater solution (2 kg of NaCl per 5 L of water). These traps were deployed at sampling sites following the CFIA’s Survey Protocol: Invasive Alien Species Forestry Trapping (2019) during the summer of 2019 and 2020. Specimens caught in the traps were removed from the solution and sent to the CFIA’s entomology laboratory. The saturated salt trap solutions were transferred to new Whirl-Pak bags and shipped to the Hanner Laboratory at the University of Guelph, Ontario, Canada. They were temporarily stored at -80˚C until further processing in the laboratory. Unfortunately, the collection fluid from two traps (trap 2 and trap 5) placed on August 26th, 2019, in Brampton could not be successfully collected. Additionally, the sample collection bag for trap number 3, collected on August 16th, 2020, in Hamilton, was empty as the salt trap solu- tion had completely evaporated at the sampling site. Insects collected during the 2019–2020 collections were not retained, while those from the 2022 collection were sent to the Hanner lab for morphological validation. Filtration The samples were taken out from the -80˚C freezer and transferred to decontaminated bea- kers, allowing them to thaw at room temperature overnight. All filtration equipment was steril- ized using a 50% bleach solution or ELIMINase from Decon Labs, Inc., followed by thoroughly rinsing with deionized water. Before filtration, the samples were inverted to resus-

pend any settled particles within the bag. Filtration of salt solutions was carried out using Nitrocellulose Mixed Ester membrane filters (47 mm diameter, 1 μm pore size) from Sterli- tech. The filtration setup consisted of a three-piece manifold connected to a GAST vacuum pump from GAST Manufactured, Inc, utilizing magnetic filtration cups from Pall. The filter cups were covered with Kimwipes throughout the procedure to prevent DNA contamination. Some samples required multiple membranes due to the presence of a high amount of debris, leading to filter clogging. The filtration time for each sample and the number of filters used PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 4 / 24

were carefully documented. As a negative control, one sample of a saturated salt solution pre- pared in the laboratory underwent filtration as well. The filters were then placed in sterile Petri dishes and stored at -80˚C until further processing for eDNA extraction. eDNA extraction Following the described steps, we performed eDNA extraction on the total number of filter membranes resulting from filtration, including negative controls. The DNeasy Blood and Tis- sue Kit from Qiagen was used with modifications to the manufacturer’s protocol [24]. Each fil- ter was thawed and divided into quarters, then further cut into strips using sterile razor blades. Small strips from each quarter were placed in individual 2 mL microcentrifuge tubes along with approximately 250 mg of glass beads (0.75–1 mm diameter). Buffer ATL (380 μL) was added to each tube, and the tubes were subjected to disruption in a Qiagen TissueLyser II for 1 minute at 30 Hz. After spinning down the samples, Proteinase K (20 μL) was added to each tube, followed by vortexing and overnight incubation at 56˚C at 700 rpm on the orbital shaker. The tubes were vortexed again, spun down, and Buffer AL (400 μL) was added. After vortexing and incubating at 56˚C for 10 minutes, ethanol (96–100%) (400 μL) was added and vortexed for another 30 seconds. All samples, including any residues, were transferred to DNeasy Mini spin columns and centrifuged at 11,000 g for one minute. The columns were then washed twice using AW1 and AW2 buffers, following the manufacturer’s protocol. After each wash,

extracts and measuring ng/uL of DNA three times (ten minutes apart) for each sample. Qubit quantifications provide a fluorescent-based estimation of DNA concentration while comparing the intensity of fluorescent dye bound to target dsDNA, in this case, to dye binding standards (DNA probes) of known concentration. As there is a dynamic in the dye intercalation and, conse- quently, variation over time in the physical interaction between the dye and the target dsDNA, estimations of DNA concentrations were taken three times and ten minutes apart. Values from each estimation were averaged, and standard deviations were calculated to increase DNA concen- tration estimation accuracy and evaluate dispersion from the mean values. A 5 μL subsample from each extract was also visually assessed for eDNA presence and quality using 1% agarose gel electrophoresis. All eDNA extracts were stored at either 4˚C for immediate use or -20˚C for short-term storage before library preparation for High-Throughput Sequencing (HTS). PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 5 / 24

Library preparation The COI gene high-throughput sequencing library was prepared through a two-step PCR pro- cess/following the same Hanner Laboratory standardized protocols and already published [4, 25]. In the first PCR, a ~407 base-pair segment of the mitochondrial genome (mtDNA) was amplified as a single fragment for all samples. The forward primers mLepF1_MiSeq/RonM- WASPdeg_MiSeq and reverse primers LepR1_MiSeq/HCO2198_MiSeq were used, along with Illumina adaptors. Each reaction included a 25 μL final volume with 2.5 μL of DNA template, 12.5 μL of 2X KAPA HiFi HotStart Ready Mix, and 0.2 μM of each primer. The cycling condi- tions consisted of 94˚C (2 min), followed by 5 cycles of 94˚C (40 s), 45˚C (40 s), 72˚C (1 min), then 35 cycles of 94˚C (40 s), 51˚C (40 s), 72˚C (1 min), and a final extension at 72˚C (5 min) [4]. The PCR products were then analyzed on a 1% agarose gel and purified using 0.8-1x NGS magnetic beads (Machery-Nagel) ratio according to the manufacturer’s instructions. For the second PCR, de novo synthesized index primers equivalent to the Nextera XT Index Kit were incorporated into the COI-Illumina adaptor sequence amplicons. This PCR step introduced unique index primer combinations for each sample. Reactions were set up in a final volume of 50 μL, containing 25 μL of 2X KAPA HiFi HotStart Ready Mix, 10 μL of molec- ular biology grade water, 5 μL of each index primer (10 μM), and 5 μL of the first cleaned-up PCR product. The PCR profile included an initial denaturation at 95˚C for 3 minutes, followed

by 8 cycles of 95˚C (30 s), 55˚C (30 s), 72˚C (30 s), and a final extension of 72˚C (5 min) [4, 26]. The PCR products were visualized on a 1% agarose gel, purified with a 0.6-1x NGS magnetic beads (Machery-Nagel) ratio, following the manufacturer’s instructions. The second PCRs were conducted in separate rooms and PCR workstations from the first PCR to minimize contamination. DNA extracts and post-PCR products were treated and pro- cessed in different areas, following standard protocols for eDNA library preparation in the Hanner laboratory. High-Throughput sequencing The Genomics Facility of the Advanced Analysis Centre (AAC) at the University of Guelph conducted the sequencing process. To prepare the sequencing libraries, the SequalPrep Nor- malization Kit from Thermo Fisher Scientific was used for normalization. The libraries were then pooled and quantified with the Qubit dsDNA HS assay kit (Thermo Fisher Scientific), while fragment size was checked using a Bioanalyzer HS DNA Chip from Agilent. Only librar- ies that passed the quality control were selected for sequencing. An Illumina MiSeq System with a MiSeq reagent kit v3 (600 cycles) was employed for sequencing, with each sample lim- ited to 1% of the run (with a maximum of 100 samples per MiSeq run). The MiSeq Reporter software was utilized for demultiplexing and adapter trimming, resulting in two paired-end raw FASTQ files per sample. Data analysis The quality control process for each raw data file (FASTQ files) involved assessing base

sequence quality scores, per base sequence content, sequence length distribution, sequence duplication levels, overrepresented sequences, and adapter content using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). Summary reports for each dataset were gen- erated using MULTIQC [27]. For COI data analysis, mBRAVE (multiplex barcode research and visualization environment) (http://www.mbrave.net/) was utilized. mBRAVE is a multi- user platform that supports the storage, validation, and visualization of High-Throughput Sequencing (HTS) data. The specific parameters set within the mBRAVE platform were the same already established in previous analyses conducted by our laboratory and proven to be PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 6 / 24

effective under similar experimental conditions [4, 25]: (1) trimming (30bp from the front and end, with a resulting trim length of 450bp); (2) filtering (minimum quality value (QV) of 20qv, minimum length of 200bp, maximum bases with low QV [<20] limited to 25%, and maximum bases with ultra-low QV [<10] limited to 5%); (3) pre-clustering threshold: none, ID distance threshold: 3%, minimum OTU size: 5, and OTU threshold: 2%; (4) assembler minimum over- lap of 20bp and maximum substitutions of 5bp; (5) The resulting sequences from this protocol were compared to libraries assembled from BOLD (Barcode of Life Data System [5]) records and matched to a BIN (alphanumeric code representing a molecular operational taxonomic unit) when possible. Each BIN represents a putative species and includes accompanying taxo- nomic information for the included records when available [28]. In addition to mBRAVE, MetaWorks was used as an alternative analysis to investigate the presence of regulated species [19]. MetaWorks is a pipeline for analyzing High-Throughput Sequencing data and is beneficial for its ability to analyze different markers and resulting sta- tistical bootstrap confidence values. MetaWorks takes raw paired-end reads as input, pair reads, then trims primers using CUTADAPT. Following trimming, MetaWorks dereplicates, and denoises reads using VSEARCH, then clusters and assigns taxonomy. An RDP (ribosomal database classifier) v2.13 was used for taxonomic classification (https://github.com/ terrimporter/CO1Classifier), accessed April 14, 2022. The parameters for MetaWorks are as

follows: (1) pairing (min quality: 20; min overlap; 25, the max fraction of mismatches: 0.02; min fraction of overlap: 0.90); (2) trimming and filtering (min length: 200bp; error rate: 0.1; end quality score: 20,20; min adapter overlap: 3; max number of N’s: (3); primer trimming for previously used COI primer cocktails reported in Milia´n-Garcı´a et al., 2020 [4]. MetaWorks outputs a results table including Exact Sequence Variants (ESVs), sequences, taxonomic assignments, and bootstrap support values for each taxonomic rank (S2 Table). In order to determine if any regulated species were present in the samples, the outputs from mBRAVE and MetaWorks were analyzed. The species returned by mBRAVE and MetaWorks (equal to and above a bootstrap value of 0.97) were compared against the updated regulated species lists in RStudio (version 022.02.0+443). To visualize the spread of regulated species, interactive and static maps were created using the leaflet, mapview, html widgets, and webshot packages in R. As metabarcoding reference databases often do not include all species of interest, the amount of missing data in the reference databases used by mBRAVE and MetaWorks was assessed. The regulated species list was compared against the databases used by both pipelines. Since mBRAVE uses BOLD, genus-specific data was accessed from BOLD on May 31, 2022, in R. Each genus from the regulated species list was mined from BOLD, which allowed for the regulated species to be compared to each species from BOLD. Since MetaWorks uses an RDP

classifier [19] (https://github.com/terrimporter/CO1Classifier), the regulated species were compared to the training set used (accessed on May 19, 2022). The effectiveness of a pipeline for IAS detection was considered based on the absolute number of molecular detections that matched morphological conformation. The higher the number of matches, the higher the effi- cacy of the bioinformatic pipeline in generating validated results. To assess potential contamination and tag-jumps, metabaR [2] was used after initial bioin- formatic analysis with MetaWorks [19]. metabaR is an R package that evaluates data quality by attempting to identify molecular artifacts using taxonomic outputs from previous bioinfor- matic pipelines, sample information, and PCR attributes. Each sampling year was analyzed separately in metabaR using the ESVs from MetaWorks and plate diagrams to determine the location of positive and negative controls. All input files were created manually, except for the MOTUs file, which was populated using a custom script, where duplicates were removed as per metabaR instructions. Data quality was assessed using each collection year’s PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 7 / 24

summary_metabarlist() and gghill_rarefaction() functions. Contaminants and tag-jumps were identified using contaslayers() and tagjumpslayer(). Additional metabaR [2] functions were utilized for data visualization. Results The COI datasets analysis of 2017 and 2018 samples generated 4,105,266 and 22,798,502 reads from 12 and 48 traps, respectively (Table 1). From the complete filter analysis of the 2019 and 2020 salt trap solutions, 14,189,292 and 39,580,980 sequence reads were generated from 23 and 46 samples, respectively (Table 1), while the quarter filter analysis of the 2020 samples yielded 45,115,854 sequence reads (Table 1), including reads observed in the negative controls (e.g., PCR negative and Sequencing controls). Similarly, 16,814,552 sequence reads were obtained for the 2022 samples, including field negatives and laboratory controls. The number of exact sequence variants (ESVs) among collection years ranged from 664 to 10,583 and com- prised 35 and 673 species, respectively (Table 1). No positive correlation was observed between the number of generated reads and the molecular operational taxonomic units (MOTUs) identified in all datasets except 2017, sug- gesting the decided sequencing depth in most cases was enough to recover the sample diversity given the primer cocktails used (Fig 1). These results were also ratified by rarefaction curves reaching the plateau for all samples in every collection year, except some in 2017 (e.g., D2; Fig 2). The sample analysis for each dataset in their experimental context showed no amplifica-

tion issues during library preparations. However, it also highlighted a considerable amount of reads in PCR controls (e.g., 2019 dataset; Fig 3), although they did not represent any MOTUs (Fig 3), suggesting spurious amplifications or PCR artifacts. No inconsistencies in PCR ampli- fication were observed for any dataset across rows or columns when analyzing the samples in their experimental context, which suggest no primer dysfunctionality or bias during library preparations. No contamination signal was observed; consequently, no contaminants were identified in the datasets based on the negative controls used. IAS detection and bioinformatics pipeline The number of IAS detected varied between mBRAVE and MetaWorks, with the mBRAVE pipeline consistently picking up a lower number of identifications in most collections except for the 2019 collection (Table 1; Fig 4). When the number of matches between molecular and morphological detection were analyzed per bioinformatic pipeline, MetaWorks resulted in the highest number of coincidences (Fig 4; Table 2). Consequently, we focused on the species that matched both bioinformatic pipelines (all except Rhagoletis cerasi and Lymantria albescens only found by mBRAVE), considering them to be “bioinformatically validated” and used MetaWorks as our reference pipeline moving forward. Table 1. The number of raw sequence reads (Seq reads) that were generated per dataset analyzed in the present study, including Exact Sequence Variants (ESVs),

number of Phyla (Phylum), Orders (Order), Genera (Genus), species (Species), and invasive alien species (IAS) detected based on eDNA after bioinformatic analysis with a bootstrap cut-off of 0.97 (MW = MetaWorks, and mBRAVE). Dataset Seq reads ESVs Phylum Order Family Genus Species IAS (MW) IAS (mBRAVE) 2017 4,105,266 664 2 12 32 35 35 1 1 2018 22,798,502 10,583 17 67 251 557 673 3 3 2019 14,189,292 5,051 13 53 204 397 445 4 5 2020 39,580,980 4,780 13 42 180 406 466 5 4 2020 Q 45,115,854 4,177 13 44 185 404 451 5 2 2022 16,814,552 3,960 9 42 182 389 444 4 3 https://doi.org/10.1371/journal.pone.0290036.t001 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 8 / 24

eDNA from up to five IAS regulated by Canada was found in the salt trap solutions with the highest bootstrap support (> = 0.97): Agrilus planipennis, Daktulosphaira vitifoliae, Lymantria dispar, Popillia japonica, and Trichoferus campestris (Table 1; Fig 4). At the time of analysis, Tri- choferus campestris was listed as a regulated pest and is therefore treated as such in the current study, however studies of recent incursions suggest this insect is not a primary pest of trees [30] and will likely be removed from the list as a result. The analysis per collection year showed that all Fig 1. Representation of MOTU numbers (nb_motus) versus the number of reads (nb_reads) per collection. Samples (NA) are represented in grey, while sequencing and PCR controls are represented in red and brown. https://doi.org/10.1371/journal.pone.0290036.g001 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 9 / 24

PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 10 / 24

five IAS together were only obtained at the sites collected in 2019 and 2020. eDNA from four of these five IAS (all but Trichoferus campestris eDNA) was detected in the 2022 sites. All IAS detected in 2022, except Lymantria dispar, were successfully identified morphologically from the decanted specimens collected in that collection year. At the sample level, there were 25 eDNA- based individual detections in 2022 (Table 2), with 13 of them lacking morphological confirma- tion. Ten out of those thirteen detections missing morphology-based conformation came from the same collection site (Windsor). The only IAS that was genetically and morphologically con- firmed to be present at Windsor was Agrilus planipennis (see sample WR_21_7_22_G3; Table 2). The remaining three mismatches occurred for Popillia japonica and Daktulosphaira vitifoliae and were only seen in two traps belonging to the same site (PB_10_8_22_3 and PB_10_8_22_4) where the species were also morphologically confirmed. Lymantria dispar and Daktulosphaira vitifoliae eDNA detections at Windsor in 2022 represented the only cases where the species could not be morphologically confirmed at any collected trap in the given site. For the 2018 collected sites, eDNA from Agrilus planipennis, Lymantria dispar, Daktulo- sphaira vitifoliae was also found, while in the samples collected from sites in 2017, only eDNA from Daktulosphaira vitifoliae was detected (Fig 4). Morphological validation In the current study, the 2022 collection was the only one where trap contents, including speci-

mens and debris, were sent to the Hanner Laboratory for morphological validation. In 2022, four IAS species were detected, with a total of 25 individual detections over three sites, 14 traps, and four collecting dates (Table 2). Of those individual detections, 11 were confirmed by both molecular and morphological methods, one was morphological only, and the remainder (13) were molecular detections only. If assessed by site and date, there is an 89% congruence between molecular and morphological data, since IAS species not detected morphologically in one trap at a site were often detected in another trap at the same site. The Windsor site collec- tions from July 21, 2022, show consistently low congruence between molecular and morpho- logical identification, with only one of 11 morphologically validated detections, which suggests a sample-specific bias. Full filter analysis vs quarters The whole filter versus quarter analysis for the 2020 collection showed that all five IAS eDNA was detected in only one quarter (Quarter I; Fig 5). These results were consistent when com- paring IAS species recovered from half and whole filters (Fig 5). Agrilus planipennis eDNA was missing in three out of the four quarters belonging to the same filters (Fig 5). Additionally, species composition analysis demonstrated that quarters or halves alone did not necessarily reflect the species composition of the whole filter, with complete filter analysis consistently showing higher average species richness for all 2020 sampled sites (Fig 6). Canopy level sampling vs ground level sampling

Only eDNA from one IAS (Lymantria dispar) was detected at the canopy level, contrary to the five IAS that were seen at the ground sampling level during the whole filter analysis of the 2020 Fig 2. Rarefaction curves per collection year (2017 and 2018) constructed with three diversity indices. D0 is species richness, D1 approaches the exponential of the Shannon index, and D2 is inverted of the Simpson index [2, 29]. Coverage refers to the Good’s coverage index [2]. Rarefaction curves per collection year (2019 and 2020) constructed with three diversity indices. 50 is species richness, D1 approaches the exponential of the Shannon index, and D2 is inverted of the Simpson index [2, 29]. Coverage refers to the Good’s coverage index [2]. Sample types (canopy and ground) are distinguished in the 2020 collection. https://doi.org/10.1371/journal.pone.0290036.g002 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 11 / 24

PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 12 / 24

collection (Fig 7). On the contrary, eDNA from all four IAS detected in 2022 were recovered from samples taken at the canopy level, and only three IAS were identified during the ground- level sampling effort that year (Fig 7). Discussion Due to their non-destructive and cost-effective nature, environmental DNA-based approaches, such as eDNA metabarcoding, have become the standard molecular method for multispecies identification. However, the routine implementation of eDNA methods in a regulatory context requires thorough validation to be broadly accepted by regulators [17]. In the present study, we pursue eDNA validation for IAS detection on CFIA-sampled salt trap collection fluids using a molecular multi-species identification approach and an experimental, morphological, and bioinformatic validation [17]. Effectiveness of eDNA metabarcoding for IAS detection eDNA from up to five invasive alien species regulated by Canada (Agrilus planipennis, Daktu- losphaira vitifoliae, Lymantria dispar, Popillia japonica, and Trichoferus campestris) was consis- tently found across five years (2017–2022) in high-risk areas sampled by the CFIA. The emerald ash borer (EAB), an invasive beetle from Asia, has caused extensive damage since its introduction to North America in 2002 [31]. The larvae of this wood-boring beetle feed on ash trees, forming tunnels under the bark and ultimately killing the trees within a few years [32]. This has led to a significant loss of ash trees, decreased biodiversity, and increased

vulnerability to other invasive plants [33–35]. The economic losses from EAB are in the bil- lions of dollars [36], making its regulation crucial [37]. In contrast, grape phylloxera, a native Fig 3. Total number of reads (nb_reads) after bioinformatic analyses and number of MOTUs (nb_motus) per collection year and sample type. NA: Represent the data generated from the samples while PCR, sequencing, and extraction indicate different control types. https://doi.org/10.1371/journal.pone.0290036.g003 Fig 4. eDNA invasive alien species detection per collection year. Barplots indicate the proportion of the IAS found per sample. A: 2017 collection; B: 2018 collection; C: 2019 collection; D: 2020 collection; and E: 2022 collection. https://doi.org/10.1371/journal.pone.0290036.g004 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 13 / 24

North American pest, exclusively feeds on grapevines and has historically caused significant vineyard losses [38, 39]. However, grafting European grapevines onto resistant rootstocks has helped manage this pest [39]. Despite being minor in Ontario, grape phylloxera remains eco- nomically significant worldwide [38, 39], emphasizing the need for early detection and popula- tion regulation. The spongy moth, introduced to Ontario in 1969 [40], defoliates various hardwood and softwood trees during outbreaks, reducing tree health and making them more susceptible to other threats. This moth is one of the most economically damaging insect pests [11], with outbreaks occurring every 7 to 10 years [40]. Spongy moths can also cause acute contact dermatitis when humans come in contact with larvae, which has implications for human health [41]. On the other hand, the Japanese beetle, native to Japan and introduced to North America in 1916 [42], then in Canada in 1939 [43], is highly destructive to turfgrass, landscape plants, and pastures. Its larvae damage grassroots, while adults heavily feed on leaves, impacting over 300 plant species in 79 families [44] and having economic impacts on over 100 of those species [42]. The Japanese beetle causes significant economic losses and requires substantial control efforts [45]. Finally, Trichoferus campestris, a velvet longhorned beetle native to Asia, was detected in Canada in 2002 [46]. T. campestris is native to Asia and feeds on healthy or mildly stressed adults in over 40 genera of broadleaf and conifer trees in its

native range; however, it prefers apple and mulberry trees [47, 48]. While it was initially Table 2. Summary of IAS detections from three sites in the 2022 collection. A = Adult, L = Larval. Specimen size is presented as an average of all lengths when multiple specimens were present. Trap Organism(s) eDNA detection Morphological detection Number of specimens Size of the specimen(s) Lifestage FE_5_8_22_4 Agrilus planipennis YES YES 1 12.2 mm A PB_10_8–22_1 Popillia japonica YES YES 3 10.8 mm A Daktulosphaira vitifoliae YES YES 2 <1 mm L PB_10_8_22_2 Popillia japonica YES YES 6 10.4 mm A Daktulosphaira vitifoliae YES YES 2 <1 mm L PB_10_8_22_4 Popillia japonica YES NO 0 n/a n/a Daktulosphaira vitifoliae YES YES 47 <1 mm L PB_10_8_22_6 Daktulosphaira vitifoliae YES YES 3 <1 mm L PB_29_8_22_1 Popillia japonica YES YES 1 11.5 mm A Daktulosphaira vitifoliae YES YES 2 < 1mm L PB_29_8_22_3 Popillia japonica YES NO 0 n/a n/a Daktulosphaira vitifoliae YES NO 0 n/a n/a PB_29_8_22_4 Daktulosphaira vitifoliae YES YES 2 <1 mm A, L Popillia japonica NO YES 1 11.2 mm A WR_21_7_22_B1 Daktulosphaira vitifoliae YES NO 0 n/a n/a WR_21_7_22_B2 Agrilus planipennis YES NO 0 n/a n/a WR_21_7_22_B3 Agrilus planipennis YES NO 0 n/a n/a Daktulosphaira vitifoliae YES NO 0 n/a n/a WR_21_7_22_G1 Agrilus planipennis YES NO 0 n/a n/a Daktulosphaira vitifoliae YES NO 0 n/a n/a Lymantria dispar YES NO 0 n/a n/a WR_21_7_22_G2 Lymantria dispar YES NO 0 n/a n/a

Daktulosphaira vitifoliae YES NO 0 n/a n/a WR_21_7_22_G3 Lymantria dispar YES NO 0 n/a n/a Agrilus planipennis YES YES 1 12.2 mm A Daktulosphaira vitifoliae YES NO 0 n/a n/a https://doi.org/10.1371/journal.pone.0290036.t002 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 14 / 24

Fig 5. Invasive alien species (IAS) detection per 2020 collected sample. Quarter filters (I, II, III, and IV), halves (I, II and II, IV), and whole filter (full) analyses are reflected, with each colour representing one of the IAS detected based on their shed eDNA in the salt trap solution. X-axis represents individual samples and the Y-axis is the proportion of IAS found per sample. https://doi.org/10.1371/journal.pone.0290036.g005 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 15 / 24

considered a primary pest of trees, recent studies suggest otherwise, indicating that it does not cause significant damage [30]. The above species detections highlight the utility of molecular biosurveillance to detect and prevent further spread. As in previous cases, an early warning for prompt actions can pro- foundly impact preserving biodiversity, significantly reducing economic losses and contribut- ing to Canadians’ well-being. Filter quarters vs. whole filter analysis The lack of consistent IAS eDNA detection on most filter quarters from the 2020 collection suggests a patchy distribution of eDNA on the filters. This can impact species detectability and constitute a sampling bias when using filter fragments alone for eDNA detection. Therefore, pooling extracts from all filter quarters belonging to the same sample before library prepara- tion [4, 49] or whole filter analysis is suggested for more accurate IAS detection and the most comprehensive inference of species diversity in a sampling site. Multiple biological replicates using whole filters is desirable in contrast to using filter fragments as a backup plan or for use as voucher samples for future analysis. Present species composition analysis beyond IAS rati- fied that quarters or halves alone did not necessarily reflect the species diversity of the whole filter. The latter might directly affect any molecular-based inference of species presence based on filter quarters, which can be common in the eDNA field [50–52].

Fig 6. Average species richness per 2020 sample collected, including whole and quarter filter analysis. QS and SA refer to specific 2020 sampling sites. Whole and half in silico refer to the ESV bioinformatic combinations for the entire filter and their halves from filter quarters. https://doi.org/10.1371/journal.pone.0290036.g006 PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 16 / 24

PLOS ONE Unveiling invasive insect threats to plant biodiversity PLOS ONE | https://doi.org/10.1371/journal.pone.0290036 August 11, 2023 17 / 24

It is important to consider that library preparation and High-Throughput Sequencing of a given sample is a sampling process itself [2]. Consequently, sequencing results mainly reflect the diversity in the processed sample under the experimental conditions rather than the site diversity. To get a more accurate representation of site diversity, the sampling effort must be comprehensive enough, there should be no primer biases or taxonomic blind spots, and all prospective species occupying the study area should be properly documented in the reference DNA databases. Failure to meet the mentioned conditions can hinder the capacity for effective eDNA IAS identification in addition to the sampling biases stated above. Experimental, bioinformatic and morphological/ecological validation The five IAS referred to above and detected across years (2017–2022) were successfully con- firmed by two independent bioinformatic pipelines (MetaWokrs and mBRAVE), increasing the accuracy of eDNA-based detections and reducing the probability of false identifications. In addition, all IAS detected in 2022, except Lymantria dispar, were successfully identified mor- phologically among the specimens sampled in that collection year, confirming species presence in the sampling area. However, not every eDNA detection at a given sample was validated with morphological confirmation. Several factors including insect biomass, time in trap fluid, and environmental conditions, can affect the sensitivity of eDNA detection methods. The amount