Lake Paringa is a moderately sized (max. depth 58 m, area 460 ha), oligotrophic lake located on the West Coast of the South Island of Aotearoa-New Zealand (43°43′10″S, 169°24′8″E; Fig. 1). The catchment (64 km²) is mainly temperate rainforest which has been largely undisturbed by human activities (Howarth et al. 2012). The West Coast is the wettest region in the country, due to the predominant westerly airflow and orographic effect of the Southern Alps. The area around Lake Paringa receives on average between 4,000–6,000 mm of rain per year (Macara 2016). The lake is culturally significant for Māori (indigenous people of Aotearoa-New Zealand). Seasonal kāinga nohoanga (small settlements) were established near Lake Paringa around 1300 CE, and the lake provided an abundant supply of food. Lake Paringa was, and still is, a noted tuna (eel) fishery, significant spawning ground and nursery for a variety of fish species, and a significant breeding area for birds (Ngāi Tahu Claims Settlement Act 1998). European access to this lake likely began around 1875 with the construction of the Haast to Paringa Cattle Track which allowed famers to move cattle along this rugged region of the West Coast. The track originally passed directly along the Eastern side of the lake following the route of the current main highway (Fig. 1). In 1965, the Haast Pass Highway was constructed, opening up this region of the West Coast for motorised travel from Central Otago to the West Coast.

Figure 1.

The Location of Lake Paringa on the West Coast of Aotearoa-New Zealand. The black dot represents the core location. Colours represent catchment land cover classes outlined in the legend. The yellow line shows the location of the Haast Pass Highway, parts of which were once the Haast to Paringa Cattle Track. Light blue lines represent streams and rivers within the catchment. LCDB5 = land cover database version 5 - https://lris.scinfo.org.nz/layer/48423-lcdb-v41-deprecated-land-cover-database-version-41-mainland-new-zealand/.

The sedimentology of Lake Paringa has been well characterised by previous studies (Howarth et al. 2012, 2021). Sediment cores from the lake are characterised by four lithofacies that form a repeating sequence related to the earthquake cycle of the Alpine Fault. The sequence includes: (i) rapidly deposited layers inferred to be co-seismic debrites and overlying turbidites formed by shaking-induced subaqueous mass wasting; (ii) stacks of rapidly deposited layers interpreted as hyperpycnites formed during periods of elevated fluvial sediment flux in response to earthquake-induced land sliding; and (iii) layered silts formed between earthquakes when the catchment was relatively geomorphically quiescent. In the present study we targeted layered silts that have accumulated post 1800 CE which are uninfluenced by deposition associated with the Alpine Fault seismic cycle.

A 6-meter-long sediment core was taken in the deepest point of the Windbag basin of Lake Paringa (58 m; Fig. 1) using a Mackereth corer (Mackereth 1958). The sediment core was transported chilled in a refrigerated trailer (4 °C) to the laboratory and stored (4 °C) for 30 days before being frozen (-20 °C). In the laboratory, the core was thawed and split in half using an electric saw to cut the plastic core tube and a sterile metal blade was used to slice the sediment horizontally across the layers to prevent vertical contamination. The core halves were photographed using a camera (Suppl. material 1: fig. S1) and logged in detail. Subsamples of sediments for DNA analysis (n = 42; Suppl. material 1: table S1) were taken from half of the core. By targeting layered silts that were deposited from more autochthonous sources, we optimised the subsampling to target sedDNA from aquatic organisms.

The top 2–3 mm of the sediment core surface were scraped off using sterile spatulas to remove potentially contaminated sediment from the splitting process. Samples for sedDNA analysis were collected by sampling the inner sediments (ranging from 4.5–9.5 g) using sterile spatulas. Subsamples were placed in 50-mL tubes and frozen (-20 °C) until DNA extraction.

Detailed visual logs and linescan imagery of the sediment core were used to correlate the new core with the original core from Howarth et al. (2012) at sub-centimetre scale using diagnostic horizons that allow sub-centimetre-scale correlations (Frith et al. 2018; Wang et al. 2020). The original master-core chronology was based on accelerator mass spectrometry measurements of the radiocarbon (14C) content of 22 terrestrial macrofossils (Howarth et al. 2012). Howarth et al. (2012) derived a calibrated calendar age for each macrofossil in OxCal 4.1 using the P_sequence depositional model.

DNA was extracted from all sediment sub-samples (n = 42) using the Lakes ABPS extraction protocol as described in Thomson-Laing et al. (2022). Briefly, sediment samples (~3–5 g) were mixed with sodium hydroxide (6 mL, 0.33M) and Tris-EDTA (3 mL, pH 8), vortexed (1 min), incubated at 65 °C for 50 min and centrifuged (15,000 × g, 1 hr, 15 °C). The resultant supernatant (7.5 mL) was transferred to a new tube and neutralized with 7.5 mL of Tris HCl (1M, pH 6.7). Sodium acetate (1.5 mL, 3M, pH 5.2) and molecular grade 100% ethanol (30 mL) were added, and samples were frozen (-20 °C) overnight to precipitate DNA. Samples were centrifuged (10,000 × g, 1 hr, 15 °C) and the supernatant discarded. DNA was extracted from the entire pellet using multiple DNeasy PowerSoil DNA Isolation Kits (Qiagen, Germany; one per 0.5 g pellet) and the automated Qiacube, following manufacturer’s protocol. Final elution volume per extracted sample was 100 µL. Eluted DNA from the same pellet extraction were pooled together following extraction. Two negative extraction controls, where no sediment was added to the extraction tube, were included every 23^rd extraction.

PCR was undertaken using primers targeting the mitochondrial 12S rRNA gene (~100 base pairs [bp]) for fish, the 18S rRNA gene (~440 bp) targeting all eukaryotes, and two regions of the mitochondrial cytochrome c oxidase I (COI) gene (~178 bp [short COI] and ~316 bp [long COI]) targeting all eukaryotes. All primers included Illumina^TM adapter tails for dual-index sequencing (Kozich et al. 2013). The 18S rRNA and 12S rRNA primers are previously published, whereas the long COI forward and short COI reverse primers were modified slightly for better recovery of zooplankton species found in Aotearoa-New Zealand (Table 1). A simple in-silico analysis was undertaken using 99 freshwater zooplankton species found in Aotearoa-New Zealand (14 cladocera, 35 copepod and 50 rotifer species) retrieved from the BOLD database (Ratnasingham and Hebert 2007). Alignment of zooplankton sequences alongside the freshwater macroinvertebrate long COI forward primer (Elbrecht and Leese 2017) and short COI reverse primer (Vamos et al. 2017) suggested the adapted primers allowed for 34 additional bp matches, resulting in only 15 bp mis-matches (EukCOI-L, BF1-adapted) and 42 additional bp matches, resulting in only 11 bp mis-matches (EUKCOI-S, fwhR1-adapted) using the newly adapted primers.

PCR reactions (25 µL) contained 12.5 µL MyFi 2 × PCR mastermix (Bioline Reagents, London, UK), 2–5 µL of template DNA and nuclease-free water, 1 µg bovine serum albumin (UltraPure BSA, Invitrogen, Massachusetts, United States) and 450 nM of forward and reverse primer. For all PCRs, temperature cycling conditions involved an initial denaturation step of 95 °C (5 min), then assay-specific cycles of denaturation, annealing and extension (Suppl. material 1: table S2) finishing with a final extension step at 72 °C (7 min). A positive control (DNA from Anguilla australis for the 12S rRNA primers and from Sabella spallanzanii for the 18S rRNA and COI primers) and a negative template control (nuclease-free water) were included in each PCR run.

To confirm successful amplification, PCR products were visualized on a 1.5% agarose gel with Nucleic Acid Staining Solution (iNtROM Biotechnology, South Korea). PCR products (20 µL) including negative extraction and PCR controls were cleaned and normalized (to ~1 ng µL^-1) with SequelPrep Normalisation Plates (Applied Biosystems, CA, USA). Cleaned PCR products were then sent to Azenta (Suzhou, China) or Sequench Ltd. (Nelson, New Zealand) for sequencing (paired-end) on an Illumina Miseq^TM with the V2 2×250 bp cycle kit for the 18S rRNA and COI amplicons and the V2 2×150 bp cycle kit for 12S rRNA fish amplicons. A bioanalyzer and Qubit 4 fluorometer (Thermo Fisher Scientific, Massachusetts, USA) were used to quantify the concentration and quality of the pooled library before it was denatured and diluted to a loading concentration of 6 pM with a 15% PhiX spike. The raw sequences were deposited into the NCBI short read archive (12S rRNA: PRJNA1007422; COI: PRJNA1007424; 18S rRNA: PRJNA1007437).

Processing of sequences for all genes was undertaken in the same way unless otherwise stated and scripts can be found at https://github.com/jkpearmanbioinf/sedimentaryDNA-trout. Samples were automatically demultiplexed on the MiSeq machine based on the dual indexing. Primers were removed from the resulting sequences using Cutadapt (Martin 2011), allowing for a single mismatch. Construction of Amplicon Sequence Variants (ASVs) was undertaken in the DADA2 package v1.26 (Callahan et al. 2016) within R version 6.2 (R Core Team 2021). Forward reads were trimmed to a length of 120 bp, 230 bp, 130 bp and 165 bp for the 12S rRNA, 18S rRNA, fwhR-amplified COI and BF1-amplified COI genes, respectively, while the reverse reads were trimmed to 120 bp, 228 bp, 130 bp and 165 bp for the 12S rRNA, 18S rRNA, fwhR-amplified COI and BF1-amplified COI genes, respectively. Sequences for each gene were denoised by a pseudo-pooling approach using an error profile constructed from a parametric model derived from the first 1×10⁸ bp of the sequences. ASVs were constructed by merging forward and reverse reads allowing for one mismatch and having a minimum overlap of 10 bp. Chimeric sequences were removed using the removeBimeraDenovo script within the DADA2 package. ASV tables are included as Suppl. material 2: tables S3–S11. Taxonomy was assigned in a two-stage process. Firstly, ASVs were classified using the RDP Naïve Bayesian Classifier (Wang et al. 2007) within DADA2 with a bootstrap of 70 for the 12S rRNA and 18S rRNA and 50 for the COI to enable classifications at higher taxonomic levels. The PR² (Guillou et al. 2013) database was used for the classification of 18S rRNA ASVs while databases were constructed based on NCBI (Agarwala et al. 2018) for 12S rRNA and BOLD (Ratnasingham and Hebert 2007) for COI made using CRABS v0.1.7 (Jeunen et al 2023). Secondly, species level classifications for the 12S rRNA dataset were determined using blastn (> 99% similarity for species classification) against the non-redundant NCBI nucleotide database. The ASVs, taxonomy and accompanying metadata were combined into ‘phyloseq’ v1.44 (McMurdie and Holmes 2013) objects which were the basis for statistical analysis.

Core sub-samples for pigment analysis (n = 20) were thawed at 4 °C, weighed in Falcon tubes and extracted three-times using acetone and a bath sonicator containing ice (30 min). The extract was dried under a stream of nitrogen gas at 40 °C and stored at -20 °C and protected from the light until analysis (< 1 week). The dried extract was re-suspended in acetone on the day of analysis and transferred to a septum-capped amber vial. Extracts were analysed by HPLC with diode array detection (DAD) using an Agilent 1260 HPLC-DAD system (Santa Clara, CA, USA) as described in Short et al. (2022) but using a shorter gradient. Extract components (10 μL injection) were loaded on to a Develosil RP Aqueous column (Nomura Chemicals; 250×4.6 mm; 5-μm) at 10%B where the solvent composition was held for 2 min before proceeding in a linear gradient to 40%B over 28 min. The column was flushed with 80%B for 3 min and re-equilibrated with 10%B for 5 min. The results from all pigments analyzed (alloxanthin, canthaxanthin, diatoxanthin, lutein, myxoxanthophyll and zeaxanthin) were summed to provide an overview of changes in total algal biomass. This approach was used, rather than using chlorophyll-a, as most of the organic matter in Lake Paringa is terrestrial (Frith et al. 2018), which makes it likely that most of the chlorophyll-a is also terrestrial and therefore not necessarily the best proxy for autogenic productivity. Mass accumulation through the core was derived from sediment bulk density and water content, where bulk density was calculated from CT densitometry using the relationship of Reilly et al. (2017). The pigment data was converted to a flux using CT-derived mass accumulation rates to account for changes in sedimentation rates through time. The pigment data was assigned as pre and post S. trutta introduction using the date established from sedDNA. To determine if there was a significant difference before and after the arrival of S. trutta, a t-test was performed after assessing variance equality on the total pigment concentration and the total pigment flux data.

All statistical analyses were performed in R (R Core Team 2021). Sequencing reads were pre-processed as follows. Sequences found in any control sample for each gene were subtracted from the datasets. Taxa with less than five total reads per dataset were removed. Sequencing datasets were rarefied to a suitable sequence depth for each gene after considering rarefaction curves (Suppl. material 1: fig. S2). The 12S rRNA gene reads were rarefied to 5,000 reads and then subset to remove birds, mammals and other taxa keeping only the Actinopteri class focusing only on fish (Fish12S rRNA). Terrestrial plant (terrestrial Streptophyta) reads were removed from the 18S rRNA gene as they were not relevant to this study, before the remaining 18S rRNA gene reads were rarefied to 20,000 reads. This gene was then subset into three groups: total eukaryotes (Euk18S), from which Chordate reads were removed, resulting in eukaryotes without fish; a subset of only fish (Fish18S), where reads were subset to include only the Actinopteri class; and groups heavily predated by trout (Prey18S) in which reads were subset to include Arthropoda, Mollusca and Rotifera. Diet studies have shown that S. trutta predate on both aquatic and terrestrial arthropods (Rutledge 1991; Sagar and Glova 1995), therefore arthropods from both habitats were included in this group. The short COI gene sequences (EukCOI-S) were rarefied to 25,000 reads and bacteria, land plants (terrestrial Streptophyta), Chordates and terrestrial worms (Megascolecidae) removed. EukCOI-S was then subset into both trout diet (Phyla: Arthropoda, Mollusca and Rotifera; PreyCOI-S) and Algae (AlgaeCOI-S). The long COI gene sequences (EukCOI-L) were rarefied to 4,000 reads and subset in the same way as the EukCOI-S, resulting in PreyCOI-L and AlgaeCOI-L.

Downcore stratigraphic plots were created for each taxonomic group analysed from the different genes (Fish12S, Fish18S, Euk18S, EukCOI-S, EukCOI-L, Prey18S, PreyCOI-S, PreyCOI-L, Algae18S, AlgaeCOI-S and AlgaeCOI-L) using the strat.plot function in the ‘rioja’ package (Juggins 2022). Constrained hierarchical clustering analysis was used to determine taxa assemblage zones downcore, constrained by sampling order to determine community shifts leading to different assemblages. A Bray-Curtis dissimilarity matrix was calculated on the plotted data using the ‘vegan’ package v2.6-4 (Oksanen et al. 2020). Hierarchical clustering was then performed using the Constrained Incremental Sum of Squares (CONISS) approach (Grimm 1987) using the ‘chclust’ function in the ’rioja’ package (v1.0-5). The primary major assemblage shift was determined; however, if these shifts occurred post-1950 they were likely associated with anthropogenic pressures such as increased human access and activity, particularly the construction of the Haast Pass Highway. In these cases, the second most significant assemblage shift was also reported. For visual purposes, taxa that showed marked changes either side of the assemblage zone were included in the plot.

Principal response curves were calculated for each taxonomic group based on Bray Curtis dissimilarity matrices. The function ‘prcurve’ from the package ‘analogue’ v0.17-6 (Simpson and Oksanen 2021) was used to calculate principal response curves for each depth downcore, using correspondence analysis to initialize the curve. Principal response curve scores were then plotted against Age and a Generalised Additive Mixed Model (GAMM) was fitted with an autocorrelation structure of order 1 (CorCAR1) to account for temporal autocorrelation using the ‘mgcv’ package v1.9-0 (Wood 2011). Due to the non-linearity of GAMMs, the method of finite differences was used to estimate derivatives of the spline to determine significant points of change in the model. Confidence intervals for the first derivative of the function were then calculated and plotted using the ‘gratia’ package (Simpson 2023), allowing for significant periods of change to be identified.

Using the 12S rRNA gene (Fish12S), nine fish taxa were found throughout the sediment core with fish reads averaging 313 reads per sample (± 113 SE). Seven of these taxa were identified to species level, six of which were native New Zealand fish (Galaxias argenteus, Galaxias fasciatus, Retropinna retropinna, Anguilla australis, Anguilla dieffenbachii and Gobiomorphus breviceps), while Galaxias sp. was only identifiable to genus level (Fig. 2a). Salmo trutta (brown trout) and Oncorhynchus tshawytscha (Chinook salmon) were also identified, appearing in 1933 (1923–1943 CE with 95% high probability density function range) and 1964 (1961–1967 CE) based on this gene, respectively.

Figure 2.

Heatmaps of fish taxa found in the Lake Paringa sediment core before and after the introduction of Salmo trutta (brown trout). Amplicon sequence variants (ASVs) taxonomically assigned to species level using the 12S rRNA primers are shown in (a), and individual ASVs from the 18S rRNA primers are shown in (b).

Using the Fish18S dataset, four ASVs from the genus Galaxias and one ASV from the Salmo genus were recovered, while five other ASVs were only classified to family (Gobiiformes) or class (Actinopteri). In this dataset, an unclassified species of Galaxias was the most relatively abundant ASV throughout the sediment core, recovered in every sample at high read abundances (Fig. 2b). Given the only Salmo found in Lake Paringa is S. trutta, we used the earliest 18S rRNA detection as the first sedDNA confirmation of this species in the lake (1906; 1892–1919 CE). Detection of S. trutta using the 18S rRNA gene was relatively consistent thereafter (Fig. 2b).

Differences in the native fish community were evident pre and post S. trutta detection, with the Fish12S data revealing Galaxias argenteus (giant kōkopu) as well as G. fasciatus (banded kōkopu) appearing only during the pre-trout phase (Fig. 2a). Salmon DNA (Oncorhynchus sp.) was only recovered in one sample from 1964 (1961–1967 CE), while S. trutta was recovered in five samples from 1933 (1923–1943 CE) to 1954 (1950–1958 CE).

After rarefaction and removal of terrestrial taxa there were on average 11,780 reads per sample and a total of 3,123 ASVs for the Euk18S dataset. The EukCOI-L dataset had an average of 10,847 reads and 3,948 ASVs, while the EukCOI-S dataset had an average of 3,813 reads and 7,521 ASVs. The Euk18S community was dominated by Fungi and Apicomplexa, with a large proportion of unclassified taxa (Suppl. material 1: fig. S3a). Both EukCOI datasets were dominated by Annelida, Arthropoda and Chlorophyta, while the EukCOI-L community also had high proportions of Fungi (Suppl. material 1: fig. S3b, c). The CONISS analysis of the eukaryotic communities retrieved using all three primer sets indicates the main shift in community structure occurred between 1863 (1852–1881 CE) and 1876 (1857–1884 CE; Fig. 3). In the Euk18S dataset, this shift was driven largely by unclassified eukaryotic taxa, an increase in chytrid fungi and shifts within the alveolate (Leidyana sp.) community. In the COI datasets (EukCOI-L and EukCOI-S), this zonation is related to shifts in unclassified eukaryotes and green algae (unclassified chlorophytes and Mychonastes sp.). GAMM analysis of both the Euk18S and EukCOI-S datasets corroborated the CONISS, with significant periods of change encompassing the CONISS-indicated shifts found for each gene (1876 [1857–1884 CE] and 1863 [1852–1881 CE] respectively; Suppl. material 1: fig. S4). However, GAMM analysis of the Euk18S dataset suggests that eukaryotic community composition was undergoing significant change over a much broader time period from 1810 (1793–1826 CE) to 1920 (1908–1933 CE; Suppl. material 1: fig. S4a), while no significant period of change was signalled from the GAMM analysis for the EukCOI-L (Suppl. material 1: fig. S4b, c).

Figure 3.

Statigraphic plot of eukaryotic sequences retrieved from (a) 18S rRNA (Euk18S), (b) long COI (EukCOI-L), and (c) short COI (EukCOI-S) primers across ~200 yrs of Lake Paringa sediment. CONISS derived dendrograms and assemblage change points (black dashed lines) are shown alongside selected taxa. The red dashed line indicates when Salmo trutta DNA appeared in the core (1906 [1892–1919 CE]) while the pale red box indicates the period S. trutta are expected to have been introduced (1870–1906 [1892–1919 CE]).

When subsetting the eukaryotic data into phyla likely to be consumed by trout (phyla: Arthropoda, Mollusca and Rotifera), there were on average 712 reads per sample and 29 ASVs for the Prey18S dataset. The PreyCOI-L dataset had an average of 401 reads per sample and a total of 536 ASVs while for the Prey-COI-S dataset there were an average of 567 reads per sample and a total of 869 ASVs. The majority of reads in these subsets were unclassifiable beyond order, particularly for the COI markers (51–63% unclassified) while mites (e.g., Maculobates sp., Sarcoptiformes and Trombidiformes) dominated the Prey18S dataset (Suppl. material 1: fig. S5). A Ceriodaphnia sp. appears in PreyCOI-L in the late 1960’s and increases to heavily dominate the group (~90% of reads in 1981 [1976–1987 CE]); however, this taxon was not identified using either the Prey18S or PreyCOI-S primers.

The main significant shift in these assemblages occurred between 1964 (1961–1967 CE) and 1972 (1971–1974 CE) in the Prey18S and PreyCOI-L datasets, and 1878 (1871–1895 CE). The second most important shifts occurred around 1891 (1877–1904 CE) and 1929 (1918–1940 CE) for the Prey18S and PreyCOI-L markers, respectively (Fig. 4). The more recent shifts (post-1965) appeared driven by the presence and significant increase in Ceriodaphnia sp. based on the PreyCOI-L. Shifts occurring pre-1900 appeared driven by an unclassified Spirostreptida (millipede) and unclassified Arthropods in the Prey18S and PreyCOI-S markers, respectively. The Prey18S community shift around 1891 (1877–1904 CE) was corroborated by a GAMM which found a significant period of change between around 1847 (1831–1863 CE) and 1909 (1985–1922 CE). The assemblage shift in 1929 (1918–1940 CE) using PreyCOI-L was corroborated by a GAMM (Suppl. material 1: fig. S6).

Figure 4.

Statigraphic plot of taxonomic groups commonly eaten by Salmo trutta in Aotearoa-New Zealand (Arthropoda, Mollusca and Rotifera). Sequences were retrieved from: (a) 18S rRNA (Prey18S), (b) long COI (PreyCOI-L), and (c) short COI (PreyCOI-S) primers over ~200 yrs in Lake Paringa sediment. CONISS derived dendrograms and assemblage change points (black dashed lines = primary, grey dashed lines = secondary) are shown alongside selected taxa. The red dashed line indicates when S. trutta DNA appeared in the core (1906 [1892–1919 CE]) while the pale red box indicates the period S. trutta are expected to have been introduced (1870–1906 [1892–1919 CE]).

Algae recovered from the Algae18S dataset had an average of 45 reads per sample and 40 ASVs. AlgaeCOI-L had an average of 274 reads per sample and 83 ASVs, while AlgaeCOI-S had on average 1,166 reads per sample and 142 ASVs in total. The algal communities found using all primer sets were heavily dominated by green algae (unclassified Chlorophyta, Choricystis sp., and Mychonastes sp.; Suppl. material 1: fig. S7). CONISS analysis suggested the major primary shift in algal community assemblage across the studied time period was around 1870 (1857–1884 CE), 1929 (1918–1940 CE) and 1946 (1940–1951 CE) according to the Algae18S, AlgaeCOI-L and AlgaeCOI-S markers, respectively (Fig. 5). Secondary shifts occurred in 1953 (1948–1956 CE) and 1883 (1871–1895 CE) using the AlgaeCOI-L and AlgaeCOI-S datasets, respectively.

Figure 5.

Statigraphic plot of Algal amplicon sequence variants retrieved from (a) 18S rRNA (Algae18S), (b) long COI (AlgaeCOI-L), and (c) short COI (AlgaeCOI-S) primers over ~200 yrs in Lake Paringa sediment. CONISS derived dendrograms and assemblage change points (black dashed lines = primary, grey dashed lines = secondary) are shown alongside selected taxa. The red dashed line indicates when S. trutta DNA appeared in the core (1906 [1892–1919 CE]) while the pale red box indicates the period S. trutta are expected to have been introduced (1870–1906 [1892–1919 CE]).

The shift occurring between 1870 (1857–1884 CE) and 1883 (1871–1895 CE) reflected in both the Algae18S and AlgaeCOI-S datasets was driven by a large increase in a Mychonastes sp. around 1886 (1874–1899 CE) when it increased up to 90% and 50% of the algal community, respectively (Fig. 5). Choricystis sp. (Algae18S) and Unclassified Chlorophyte (AlgaeCOI-S) ASVs also appear related to this shift. Later shifts from 1929 (1918–1940 CE) to 1946 (1940–1951 CE) captured by the AlgaeCOI-L and AlgaeCOI-S markers were driven by Sphaeropleales sp. and the decline of unclassified Chlorophyte ASVs (Fig. 5). GAMM’s applied to the Algae18S dataset suggest a significant period of change between 1883 (1871–1895 CE) and 1906 (1892–1919 CE), which given the 95% high density probability function for the age model at this time may include the CONISS delimiter around 1870 (1857–1884 CE). However, no significant periods of change were found when applying GAMMs to either of the COI datasets (Suppl. material 1: fig. S8).

The total pigment flux was highest (> 1.5 µg cm^-2 y^-1) in 1774 CE (1754–1794 CE), 1821 (1798–1844 CE), and 1970 (1968–1971 CE; Suppl. material 1: fig. S9a). The concentration of pigments was relatively low throughout the study period, with the highest concentration being 4.16 µg/g in 1774 (1754–1794 CE; Suppl. material 1: fig. S9b). A t-test showed no significant difference in total pigment flux before and after S. trutta arrival (p = 0.7, Suppl. material 1: fig. S9a).

Fish community

Salmo trutta was observed using the Fish18S marker, with sequences attributed to a Salmo sp. occurring around 1906 (1892–1919 CE), and consistent recovery thereafter. We presume that trout introduction occurred sometime after 1870 CE and speculate that biomass was insufficient for it to be detected using the molecular approach until 1906 (1892–1919 CE), as previous studies have shown that DNA detection increases with higher biomass (Nelson‐Chorney et al. 2019; Kuwae et al. 2020). Unfortunately, there are no written records to confirm the exact timing of trout introduction in Lake Paringa.

Changes in the historical fish community of Lake Paringa were evident using both 18S (Fish18S) and 12S rRNA (Fish12S) markers, with a shift from a galaxiid to S. trutta dominated community after about 1906 (1892–1919 CE). Salmonids such as S. trutta have been shown to negatively impact native fish populations across Aotearoa-New Zealand, proving particularly detrimental to galaxiid species which historically dominated many freshwater ecosystems (Rowe et al. 2002; McDowall 2006). Both markers reveal a pre-trout native fish community dominated by Galaxia sp. including G. argenteus (giant kōkopu) and G. fasciatus (banded kōkopu) which declined in presence and relative abundance just prior to the detection of S. trutta DNA (around 1906 [1892–1919 CE]) and disappear thereafter. Similar declines have been recorded in streams and rivers throughout Aotearoa-New Zealand (Crowl et al. 1992; McDowall 2006), where only one of 30 known galaxiid species is currently classified as ‘not threatened’ (Dunn et al. 2017).

Studies have shown that both predation and competition for food and habitat by introduced Salmonids, in particular S. trutta, have led to disjunct distributions and dwindling numbers of Galaxiids in streams and rivers across Aotearoa-New Zealand (Taylor and Main 1987; Minns 1990; Bonnett and Sykes 2002; David et al. 2002; McIntosh et al. 2010). In the present study, the sedDNA suggests this may also be the case for lake populations, which has been highlighted by McDowall (2006). However, the often-patchy recovery of fish sedDNA (Capo et al. 2021), such as was found for G. fasciatus, A. australis and Onchorhynchus tshawytscha using the 12S rRNA marker, indicates that more replication (e.g. cores, DNA extraction, PCR), larger sample sizes, and sediment cores taken from locations where fish are abundant (i.e. near the shore; Thomson-Laing et al. 2023) may be required to obtain greater quantities of fish sedDNA, in order to be confident results are not skewed by false negatives.

The introduction of new apex fish predators is often associated with foodweb trophic cascades through the creation and release of predation pressures at differing levels of the foodweb (Carpenter et al. 2001; Tronstad et al. 2010; Strock et al. 2013; Milardi et al. 2016). However, specific cascading effects are dependent on the number of pre-existing trophic levels occupied e.g., whether previous apex predators were piscivorous or planktivorous (Su et al. 2021). Another potential response is that in which a newly introduced predator in part or wholly replaces pre-existing predators with no resulting cascading effects. In Aotearoa-New Zealand, trout are known to displace native generalist fish taxa such as kōaro and kōkopu (McIntosh et al. 2010), therefore the lack of cascading effects on smaller, planktivorous fish (e.g., bullies, smelt and other galaxiids) may in part be due to S. trutta replacing the pre-existing generalist/piscivorous controls on these species.

Post 1906 (1892–1919 CE), a signal for S. trutta was consistently observed in the sediment core using the Fish18S marker. Salmo trutta sequences were also recovered using the Fish12S marker, although these were only found in < 20% of samples post-introduction. A comparison of these two markers shows potential similarities in fish communities; however, precise comparisons are difficult due to the disparity in taxonomic resolution between the two genes. While sequences recovered using Fish18S were only identifiable to class, family, and genus level, the Fish12S data could be assigned to species level using a manually constructed database of known 12S rRNA fish sequences. There is agreement between the markers at the class, family, and some cases genus level, for example, sequences from both markers include Galaxias and Salmo genera as well as the Gobiiformes family. The lack of taxonomic discrimination using the 18S rRNA gene could be due to this region of the gene being too conserved for fish species and/or a lack of fish sequences in the 18S rRNA database. Consequently, while this gene was successful in recovering fish DNA more consistently, particularly for S. trutta, the lack of taxonomic resolution is limiting. In the present study, the use of multiple markers proved to be complementary. We were able to capitalise on the taxonomic resolution of the 12S marker while benefiting from the more consistent recovery of fish sequences using the 18S marker. However, given that the species-specificity of metabarcoding assays can be low, another approach to identify the timing of S. trutta introduction would be the use of a quantitative, species-specific assay. A targeted S. trutta assay using qPCR or ddPCR could help improve detection probability, potentially providing a more accurate estimate of the timing of S. trutta introduction.

The cascading effects of trout introduction are often visible in aquatic ecosystems via shifts in zooplankton and algal dynamics or community structure, and even changes in zooplankton body-size distributions (Tronstad et al. 2010; Strock et al. 2013; Milardi et al. 2016). The most significant eukaryotic community shift across the ~200-year period from all markers was between 1863 (1852–1881 CE) and 1876 (1857–1888 CE; Fig. 3), encompassing the probable timing of S. trutta introduction. Interestingly, this was related to different taxa in the two different genes analysed. Using the two COI markers, green algae (Chlorophyta and Mychonastes sp.) appeared to drive the community shifts, while unclassified Eukaryotes and shifts in the relative abundances of different Alveolates (Leidyana sp.) appeared responsible for the shift detected using the 18S rRNA marker. Chytrid fungi (Chytridiomycetes) heavily dominated the Euk18S dataset, increasing in relative abundance from around 1830 (1814–1847 CE), coinciding with European settlement of Aotearoa-New Zealand. Some Chytrid fungi are highly parasitic to phytoplankton species (Kagami et al. 2007), which could explain the decline in two Choricystis sp. between around 1830 (1814–1847 CE) and 1883 (1871–1895 CE) and/or a Desmodesmus sp. which disappears in 1866 (1852–1881 CE). Due to a lack of sequences in the databases, other taxa related to this shift were largely unclassified Eukaryotes and unclassified Chlorophytes.

Lake Paringa was chosen for this study due to its relatively unmodified catchment which is almost entirely indigenous forest. While much of the country was experiencing deforestation and the rapid growth of pastoral agriculture during the time of European settlement (mid to late 1800s), the Lake Paringa catchment remained forested. European settlers began to access the surrounding area around 1875 when a small cattle track was constructed to the east of the lake, allowing small herds of cattle to be moved north. Minor amounts of nutrients may have been deposited near, or in, the lake during the movement of cattle, although the impact of these sporadic cattle-drives is likely to be low. Similarly, the introduction of wild terrestrial mammals including red deer and marsupials such as possums to the region in the mid to late 1800s may have released slightly more nutrients into the surrounding land with potential to enter the lake, but again, the in-lake effects are likely to be very minor. Algal pigment analysis corroborates the lack of productivity around that time, with only a small spike in 1850. The most major change in catchment land use would have been in 1965 when the Haast Pass Highway was constructed along the south-east side of the lake. Earthworks would have occurred and subsequent increases in vehicle traffic and visitors to the area may have impacted the lake, however this occurred much later than the introduction of S. trutta, which is the time of interest in this study. Therefore, there is a high likelihood that this significant shift in the overall eukaryotic community is related to the introduction of trout, given that no other major disturbances occur at this time that could explain the change in this broad range of taxa.

Trout diet studies reveal a broad range of organisms are consumed by these generalists including fish, aquatic and terrestrial arthropods, molluscs and rotifers (Rutledge 1991; Sagar and Glova 1995; Olsson et al. 2006). While changes in terrestrial arthropod communities may be related to terrestrial based processes such as the introduction of land mammals, they were included in this study due to the possibility that trout could reduce their sedDNA signal in the lake samples, but primarily because the majority of arthropod sequences were unclassifiable. This makes it difficult to determine which of the unclassified taxa are terrestrial, or not part of the trout diet. As such, all arthropods were included in the dataset and their habitat noted where known. It should therefore be noted that the changes in arthropod signals over time cannot be directly associated with specific trout interactions.

When exploring our data on eukaryotic groups likely to be predated by trout (Arthropods, Molluscs and Rotifers), CONISS analyses identified the primary assemblage shift in two markers occurred post-1950 (1946–1953 CE), likely related to anthropogenic changes including increased human access. Primary and secondary shifts also occurred between 1878 (1871–1895 CE) and 1891 (1877–1904 CE), for the PreyCOI-S and Prey18S dataset respectively coinciding with the probable timing of S. trutta introduction. GAMM analysis corroborated this significant period of change in the Prey18S dataset but not using either PreyCOI datasets. The lack of detection of a significant period of change in the PreyCOI subsets could be due to the inconsistent and patchy recovery of many arthropod, mollusc and rotifer taxa making temporal changes difficult to assess. Most sequences in these datasets were unable to be taxonomically classified, impeding interpretations on how trout are possibly affecting the lake community.

Different prey taxa were also found to be dominant depending on the marker used. For example, freshwater mussels (classified as Elliptio sp. but likely to be Echyridella sp.) dominated the Prey18S dataset along with a range of mites classified to genus and sometimes order level. Whereas, no specific ASVs dominated either COI datasets until 1967 (1964–1968 CE) when a Ceriodaphnia sp. appeared and increased dramatically in the PreyCOI-L dataset. The recovery of ASVs was patchy using both COI markers where many taxa were found only a few times throughout the core, making time-series analysis of these results difficult.

Of note are the lack of known Aotearoa-New Zealand zooplankton (Rotifers, Copepods and Cladocera) retrieved throughout the sediment core which are consistently lacking across all three markers except for the significant increase in a Ceriodaphnia sp. in 1967 (1964–1986 CE) using the EukCOI-L marker. This was unexpected for the two COI markers, given that in this study both the EukCOI-L forward and EukCOI-S reverse primers were specifically adapted to increase base-pair matches with a range of Aotearoa-New Zealand zooplankton species. In addition, the 18S primer has successfully retrieved zooplankton reads from contemporary water samples in a previous study (Pearman et al. 2022) indicating this may be a methodological challenge relating to DNA preservation or recovery of these organism from lake sediments. Consequently, any potential zooplankton effect could not be determined using this sedDNA approach. As zooplankton are known to be components of trout diets, the poor retrieval of these species in the DNA could potentially result in shifts in these vital components of the food web being missed and further work to improve the detection of zooplankton in lake sediment eDNA studies is required.

There was evidence of a significant shift in algal taxa around the time of trout introduction. The Algae18S marker suggests this as the most major assemblage shift, while the AlgaeCOI-S marker found this was the second most major shift. Shifts in diatom species assemblages have been recorded following trout introduction (Sienkiewicz and Gąsiorowski 2016; Cantonati et al. 2021). However, in this study, a general decrease in certain species of green algae alongside a significant increase in Mychonastes sp. and other unclassified algal species appear responsible for the assemblage shift occurring after the introduction of S. trutta. The sudden presence and subsequent increase in Mychonastes sp. found by both Algae18S and AlgaeCOI-S markers, and potentially by the AlgaeCOI-L marker as an ASV of the order Sphaeropleales, appear to have driven the significant algal assemblage shift between 1870 (1857–1884 CE; Algae18S) and 1883 (1871–1895 CE). The appearance and growth of the small (< 3 µm), coccoid, eukaryotic picoplankton Mychonastes sp. (Liu et al. 2020) could suggest a relaxation of grazing pressure at the micro- or macro-zooplankton level, following the introduction of S. trutta. However, these impacts are unable to be directly determined due to the lack of zooplankton recovered throughout the core.

The shifts in algal community structure observed in Lake Paringa did not correspond with a shift in algal pigment concentrations (as an indication of primary productivity in the lake). There were no significant changes in total pigment flux associated with this shift in algal community structure and patterns in algal pigments were consistent pre and post introduction of S. trutta. This is different from observations in other studies (Tronstad et al. 2010; Strock et al. 2013), although those algal shifts recorded as a result of trophic cascades may be due to a corresponding increase in nutrient levels in those systems. Alternatively, the lack of obvious pigment shifts could also be due to the low pre-existing algal abundance in this oligotrophic lake.

Overall, sedDNA results suggest changes to the Lake Paringa fish, eukaryotic and algal communities occurred around 1856 (1841–1871 CE) to 1891 (1877–1904 CE), which could be related to the timing of S. trutta introduction to the lake, which was first detected in the sedDNA in 1906 (1892–1919 CE). The metabarcoding approach identified possible responses from a wide range of organisms, including microscopic (e.g., fungi) and soft-bodied (e.g., polychaetes) taxa which would not be detectable using traditional paleolimnological approaches. However, there are significant caveats to the use and interpretation of sedDNA that need to be considered.

In this study, the data presented are from a single core location from one lake and are therefore limited in their potential extrapolation to other lakes. Collecting cores from multiple unimpacted lakes with introduced trout would add evidence to the pattern observed here. To explore shifts in foodweb components, we used a multiple markers approach to reduce the chance that key taxa were missed due to primer biases (Zhang et al. 2018). However, we still failed to robustly detect some taxonomic groups of interest (e.g., zooplankton). The 18S rRNA primers (Zhan et al. (2013) have been successful in previously detecting zooplankton species in Aotearoa New Zealand both in water and lake surface sediments (Pearman et al. 2022), however very few sequences were recovered in this study. Although three different markers were used, zooplankton taxa reads were absent or very low throughout the majority of the sediment core. While two of the three primers were adapted to increase recovery of zooplankton found in Aotearoa-New Zealand, further in vitro analysis needs to be undertaken to determine why these taxa were not amplified. This could be achieved through sediment-spiking, in which zooplankton remains are added to samples of core sediment before extracting DNA in order to assess DNA recovery (Brasell et al. 2022). It is possible that lake-specific sediment chemistry may negatively impact DNA recovery of certain organisms, and this should be further investigated. It should also be noted that for the algal subset of the 18S rRNA data, very low read numbers were obtained. This could be due to the universal nature of the primer targeting other eukaryotes.

The patchy recovery of certain taxa such as fish is also of concern. A number of species likely to have inhabited the lake throughout the study period were only found sporadically. For example, Anguilla australis and A. dieffenbachii, two native eels known to inhabit freshwaters of the West Coast region prior to 1700 CE, and which still exist today, are recovered in only five samples collectively. Inconsistent recovery of DNA from larger, motile organisms such as fish is a challenge in lake sediments (Thomson-Laing et al. 2022), as DNA can be heterogeneously distributed throughout lakes. However, taxa recovered using both COI primers were noticeably patchy in their recovery, meaning the community shifts calculated for those data subsets were largely driven by changes in rare and / or inconsistently recovered taxa. This could be due to the COI primers being designed to target multicellular organisms which have a more patchy distribution in sediment compared to microbial species which are also amplified with the 18S rRNA primers. Variability in DNA recovery could also be related to changes in the sedimentology, and further study is needed to explore this in detail.

While work is ongoing to populate databases with the huge influx of newly discovered sequences, these databases are still lacking in many areas, particularly those of smaller, less charismatic organisms. The overwhelming majority of ASVs recovered from this sediment core were unclassifiable below order level, illustrating the lack of Aotearoa-New Zealand species represented in major databases. While information is still contained in this data, validation and ground-truthing requires taxonomic classification, which was unobtainable for the majority of sequences. Hence the pressing need to incorporate species from diverse geographic areas such as Aotearoa-New Zealand into sequence databases.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This research was funded by the New Zealand Ministry of Business, Innovation and Employment research programme – Our lakes’ health; past, present, future (grant no. C05X1707).

Conceptualization: DK, SAW, MJV. Data curation: JKP, LAS. Formal analysis: LAS, GTL, AR, JP. Funding acquisition: MJV, SAW. Methodology: JDDH. Writing - original draft: LAS. Writing - review and editing: JDDH, GTL, JKP, DK, JP, SAW, MJV, SF, AR.

Lena A. Schallenberg https://orcid.org/0000-0003-1036-4855

Georgia Thomson-Laing https://orcid.org/0000-0001-8337-5489

John K. Pearman https://orcid.org/0000-0002-2237-9723

Jamie D. Howarth https://orcid.org/0000-0003-4365-0292

Sean Fitzsimons https://orcid.org/0000-0003-0513-7367

All of the data that support the findings of this study are available in the main text or Supplementary Information.