Kenai National Wildlife Refuge Aquatic Invasive Fish Surveys

Study Extent

Our target universe was all waterbodies in the study area susceptible to invasion by non-native fish, particularly invasive northern pike. Our initial sample frame was the set of lakes in the vicinity of the Kenai National Wildlife Refuge. We considered individual lakes to be the sampling units.

Sampling

We selected 16 lakes to survey for northern pike in 2023, mostly basing our selections on the prioritization of the Alaska Department of Fish and Game’s Invasive Species Lake Prioritization (Alaska Department of Fish and Game, 2022). We also took into account recent pike surveys, avoiding lakes that had been surveyed for pike in the last 10 years or where surveys are planned for 2024. We collaboratively planned our survey schedule with our partners. To estimate acreages of the littoral zone as defined by Dunker et al. (2022) as all of the lake area with a depth less than 4 m, we referred to available bathymetric maps. For lakes where no bathymetric maps were available, we examined aerial and satellite imagery to estimate acreage of the littoral zone for each lake. We allocated samples across the selected lakes by first allocating 5 samples to each lake, then we allocated the rest of the samples with the number of samples being proportional to estimated littoral acreages of the lakes, yielding sample sizes of 6–14 samples per lake. To decide when to sample each lake, we sorted our intended order of sampling lakes based on the pike priority values from the Alaska Department of Fish and Game’s Invasive Species Lake Prioritization, choosing to sample highest priority lakes first. In order to to minimize disturbance of nesting swans, we adjusted our schedule and selection of lakes based on results of spring swan surveys. Where nesting swans were present, we removed these lakes from our set of lakes to be surveyed and substituted other lakes where nesting swans were not present. We conducted eDNA surveys for northern pike in May–June 2023 as early in the season as was feasible for two reasons: First, northern pike are likely most detectable by eDNA methods in the spring due to spawning behavior, shedding more DNA into the water in the spring than at other times of the year (Dunker et al., 2022). Second, we wanted to send off samples early in the season so that results were available before the end of the field season, enabling us to follow up any potential positive detections with gillnet surveys. Within each lake to be surveyed, we selected sampling locations before going out into the field using Google Earth Pro (https://www.google.com/earth/versions/#earth-pro), spreading the number of sampling intended locations over the littoral zone following the guidance of Dunker et al. (2022).

Method steps

1. We collected water samples using eDNA water sampling kits provided by Jonah Ventures (Boulder, Colorado) following the kits’ sampling instructions (Jonah Ventures, 2022). At each site, using gloved hands, we first drew up 60 ml into a syringe from a depth of 1 to 15 cm. We pushed this water through a 25 mm diameter, 1 µm Entegris nylon syringe filter. We repeated these steps three times for a total of 180 ml or until the filter became clogged, resulting in water samples ranging in volume from 60 ml to 180 ml. We dried the filters by using the syringe to push air through them. We then filled the filter cartridges with Triton X-100 preservative. Samples were kept cool in a refrigerator until they could be shipped out for processing.
2. Using an R, version 4.2.3 (R Core Team, 2023) script, we randomly assigned 2/3 of the samples to be processed by a qPCR assay designed to detect northern pike (Dunker et al., 2016) and 1/3 of the samples to be processed by fish metabarcoding (Miya et al., 2015). All samples were shipped to Jonah Ventures for processing.
3. For qPCR samples, sample filters, lysis buffer, and proteinase K were heated to 56 °C for one hour. Under a laminar flow hood, warm lysis buffers were pushed through the filter housing and all supernatant was collected in the corresponding lysate tube. Tubes were placed in an incubator overnight at 56 °C. After incubation, the lysate tubes were immediately processed. Genomic DNA from samples was extracted using the DNeasy Blood & Tissue Kit (250) (catalog number 69506) according to the manufacturer’s protocol. Whole filters were used for genomic DNA extraction. Genomic DNA was eluted into 200 µl and frozen at -20 °C. An amplicon from the cytochrome oxidase, subunit I (COI) gene was amplified via qPCR from genomic DNA samples using Northern Pike COI FWD (5’ CCTTCCCC CGCATAAA TAATATAA 3’) and REV (5’ GTGTTGAA GCTGGTGC TGGTAC 3’) primers, and Northern Pike COI Probe (5’ /56-FAM/ CT+TC+TG+AC+TT+CTC+CCC/ 3IABkFQ/ 3’) of Dunker et al. (2016). A standard curve was generated for each run to correspond to targeted region of the northern pike COI gene. Each qPCR reaction was run in triplicate and contained 8.0 µl of QuantaBio PerfeCTa qPCR ToughMix Low ROX (catalog number 97065-966), 500 nM of each primer, 300 nM of probe, 4.0 µl of gDNA, and 4.8 uL of nuclease-free water for a total reaction volume of 20 µl. qPCR amplification was carried out on the Thermofisher QuantStudio 5 qPCR instrument with the following thermal profile conditions: 1 cycle of initial denaturation for 5 minutes at 95 °C followed by 50 cycles of 15 seconds at 95 °C and 1 minute at 60 °C. A standard curve was tested in triplicate for each qPCR run. The targeted gene of interest, obtained from NBCI, was used to design a synthetic gBlock (TTCCCCTA ATGATTGG TGCCCCCG ACATGGCC TTCCCCCG CATAAATA ATATAAG CTTCTGAC TTCTCCCC CCCTCCTT TTTACTTC TCTTAGCC TCCTCAGG TGTTGAAG CTGGTGCT GGTACTGG CTGAACAG TTTATCC GCCTTTGG CCGG, from Integrated DNA Technologies) that contained the northern pike primers and probe. A 10-fold dilution was carried out on the northern pike gBlock, generating a 7-point standard curve ranging from 5,372,000 to 5.372 copies. Each qPCR reaction contained 8.0 µl of QuantaBio PerfeCTa qPCR ToughMix Low ROX (catalog number 97065-966), 500 nM of each primer, 300 nM of probe, 2.0 µl of northern pike gBlock, and 6.8 µl of nuclease-free water for a total reaction volume of 20 µl. qPCR was carried out on the Thermofisher QuantStudio 5 qPCR instrument with the following thermal profile conditions: 1 cycle of initial denaturation for 5 minutes at 95 °C followed by 50 cycles of 15 seconds at 95 °C and 1 minute at 60 °C. Analysis of qPCR data was carried out using the Thermofisher Connect™ cloud software using default settings. A linear regression was applied to the calibration curve which showed the relationship between the log10-transformed standard concentration and the number of PCR cycles at which the detection threshold was reached (Cq). The R2 intercept and slope of the linear regression were examined for goodness of fit, with an R2 value >0.99 and a reaction efficiency (E), or how close to a doubling of product was achieved with each PCR cycle, within 85%–110%. A 100% efficiency is a slope of ~3.3 cycles per 10-fold dilution. Sample quantities were extrapolated from the standard curve linear regression based on the Cq value at which the detection threshold was reached and back calculated to number of copies / 100 ml in the original sample volume. More complete qPCR methods are available from Jonah Ventures (2023b).
4. For fish metabarcoding samples, a customized one-ton arbor press along with a removable leather punch was used to open the plastic casing of each filter. Once the plastic casing was cut, sample barcodes were recorded and assigned to wells within a 96 well plate or numbered extraction tubes. The whole filter was removed and transferred to the extraction plate or tube using sterilized tweezers inside a laminar flow hood. The removable leather punch was sterilized between each eDNA filter. Plates or tubes were immediately processed or stored at -20 °C until the extraction process could be performed. Genomic DNA from samples was extracted using the DNeasy Blood & Tissue Kit (250) (catalog number 69506) according to the manufacturer’s protocol. Whole 25 mm filters were used for genomic DNA extraction. Genomic DNA was eluted to 200 µl and frozen at -20 °C. Portions of hyper-variable regions of the mitochondrial 12S ribosomal RNA (rRNA) gene were PCR amplified from each genomic DNA sample using the MiFishUF (GTCGGTAA AACTCGTG CCAGC) and MiFishUR (CATAGTGG GGTATCTA ATCCCAGT TTG) primers with spacer regions (Miya et al., 2015). Both forward and reverse primers also contained a 5’ adapter sequence to allow for subsequent indexing and Illumina sequencing. PCR amplification was performed in replicates of six and all six replicates were not pooled and kept separate. Each 25 µl PCR reaction was mixed according to the Promega PCR Master Mix specifications (Promega catalog number M5133, Madison, Wisconsin) which included 12.5 µl Master Mix, 0.5 µM of each primer, 1.0 µl of gDNA, and 10.5 µl DNase/RNase-free water. DNA was PCR amplified using the following conditions: initial denaturation at 95 °C for 3 minutes followed by 45 cycles of 20 seconds at 98 °C, 30 seconds at 60 °C, and 30 seconds at 72 °C, and a final elongation at 72 °C for 10 minutes. To determine amplicon size and PCR efficiency, each reaction was visually inspected using a 2% agarose gel with 5 µl of each sample as input. Amplicons were then cleaned by incubating amplicons with Exo1/SAP for 30 minutes at 37 °C following by inactivation at 95 °C for 5 minutes and stored at -20 °C. A second round of PCR was performed to complete the sequencing library construct, appending with the final Illumina sequencing adapters and integrating a sample-specific, 12-nucleotide index sequence. The indexing PCR included Promega Master mix, 0.5 µM of each primer and 2 µl of template DNA (cleaned amplicon from the first PCR reaction) and consisted of an initial denaturation of 95 °C for 3 minutes followed by 8 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s. Final indexed amplicons from each sample were cleaned and normalized using SequalPrep Normalization Plates (Life Technologies, Carlsbad, California). 25 µl of PCR amplicon was purified and normalized using the Life Technologies SequalPrep Normalization kit (catalog number A10510-01) according to the manufacturer’s protocol. Samples were then pooled together by adding 5 µl of each normalized sample to the pool. Sample library pools were sent for sequencing on an Illumina MiSeq (San Diego, California) at the Texas A&M Agrilife Genomics and Bioinformatics Sequencing Core facility (College Station, Texas) using the v2 500-cycle kit (catalog number MS-102-2003). Necessary quality control measures were performed at the sequencing center prior to sequencing. Raw sequence data were demultiplexed using pheniqs, version 2.1.0 (Galanti Shasha and Gunsalus, 2021), enforcing strict matching of sample barcode indices (i.e, no errors). Cutadapt, version 3.4 (Martin, 2011) was then used remove gene primers from the forward and reverse reads, discarding any read pairs where one or both primers (including a 6 bp, fully degenerate prefix) were not found at the expected location (5’) with an error rate < 0.15. Read pairs were then merged using vsearch, version 2.15.2 (Rognes et al., 2016), discarding resulting sequences with a length of < 130 bp, > 210 bp, or with a maximum expected error rate > 0.5 bp (see Edgar and Flyvbjerg, 2015). For each sample, reads were then clustered using the unoise3 denoising algorithm (Edgar, 2016) as implemented in vsearch, using an alpha value of 5 and discarding unique raw sequences observed less than 8 times. Counts of the resulting exact sequence variants (ESVs) were then compiled and putative chimeras were removed using the uchime3 algorithm, as implemented in vsearch. For each final ESV, a consensus taxonomy was assigned using a custom best-hits algorithm and a reference database consisting of publicly available sequences from GenBank (Benson et al., 2005) as well as Jonah Ventures voucher sequences records. Reference database searching used an exhaustive semi-global pairwise alignment with vsearch, and match quality was quantified using a custom, query-centric approach, where the % match ignores terminal gaps in the target sequence, but not the query sequence. The consensus taxonomy was then generated using either all 100% matching reference sequences or all reference sequences within 1% of the top match, accepting the reference taxonomy for any taxonomic level with > 90% agreement across the top hits. More complete fish metabarcoding methods are available from Jonah Ventures (2023a).
5. We reshaped the qPCR and fish metabarcoding data into occurrence data suitable for publication to GBIF (https://www.gbif.org/) using a Quarto (https://quarto.org/) document (Bowser, 2023) that ran R, version 4.2.3 (R Core Team, 2023) in RStudio, version 2022.12.0. We used the R packages Biostrings, version 2.66.0 (Pagès et al., 2022); reshape2, version 1.4.4 (Wickham, 2007); and sf, version 1.0-12 (Pebesma, 2018; Pebesma and Bivand, 2023). We had used the uuid package, version 1.1-0 (Urbanek and Ts’o, 2022) to generate UUIDs before the document was rendered. Among the metabarcoding data were unexpected reads of Carassius (94 reads of one ESV from one sample and 46 reads of a second ESV from a second sample, the samples coming from different lakes) Lepomis (10 reads of one ESV from one sample), Micropterus (50 reads of one ESV from one sample), and Gadus (32 reads of one ESV from one sample). We believed that these unexpected records likely came from contamination or indexing errors, so we filtered out all occurrences with read abundances < 100 reads. We did not detect Esox lucius, our primary target species, so we added inferred absences of this species for all sampling events.

We surveyed for all non-native fish that could occur in this area, but with a particular focus on northern pike.

1. Chordata
   common name: fish rank: phylum
2. Esox lucius
   common name: northern pike rank: species

The geographic extent included freshwater lakes in the vicinity of the Kenai National Wildlife Refuge, Kenai Peninsula, Alaska, USA.