This study examines the fish biodiversity within a single catchment of the Murray-Darling Basin (MDB) (Australia). River morphology, abiotic water conditions and fish biodiversity vary along an altitudinal gradient in the MDB (Lintermans 2007). High altitude river systems are characterized by rocky substrates, high flow rates, clear waters, and low fish biodiversity. Low altitude river systems generally have a soft bottom substrate, are slow-flowing, turbid and contain higher fish biodiversity. Environmental DNA sampling was conducted from the upper section of the Murrumbidgee River (MR) catchment at five sites within the main channel. Sampling sites were selected based on their altitude and expected species diversity (Figure 1) (Lintermans 2002). Sampling sites ranged from high (1,303 m above sea-level at MR01) to low (520 m above sea-level at MR05) altitude with a minimal distance of ca. 35 km between sampling sites (measured along the river channel). The average river width at each sampling site was 11.75 m (MR01), 33.75 m (MR02), 58.25 m (MR03), 52.00 m (MR04) and 52.25 m (MR05).

Figure 1.

Map of all sampling sites within the main channel of the Murrumbidgee River (MR). Sampling sites are numbered from upstream (MR01) to downstream (MR05).

To obtain a measure of expected fish biodiversity, which is required to evaluate the overall performance of the eDNA metabarcoding survey, an expert opinion survey was conducted. Additionally, conventional boat electrofishing data was available for the three most downstream sites. A detailed description of the expert opinion survey data and the electrofishing data is available in the supplementary materials (Suppl. material 1).

Environmental DNA samples (i.e. twelve 2 L water samples) were collected from each site between the 3^rd and 8^th of November 2016 (only eleven samples were available from the MR01 site as one of the sampling bottles broke). Potential contaminating DNA was removed from sampling equipment prior to collecting water samples, using a 20% bleach solution and thoroughly rinsing with UV-sterilized tap water. One blank field control (BFC) was included for each site and consisted of a 2 L sampling bottle filled with UV-sterilized water which was opened on site, closed and submerged in the water. At each site, samples were collected over four transects across the river width spanning a 100 m river section (i.e. ca. 33 m distance between transects). At the MR01 site samples were collected by wading in the river while for the remaining sites all samples were collected from a canoe. Along each transect surface water samples were collected from both the left and right river banks (i.e. within 1m of the river bank) and the mid-channel. Samples were stored on ice, transported to the University of Canberra (ACT, Australia) and eDNA was captured within 12 h using a 1.2 µm glass fibre filter (Sartorius, Göttingen, Germany). Filtering equipment was cleaned between samples as described above and negative equipment controls (NEC) were obtained by filtering 500 mL of UV-sterilized tap water prior to processing water samples. Filters were stored at -20 °C until eDNA extractions were performed at the trace DNA laboratory (University of Canberra) using the PowerWater DNA Extraction Kit (MoBio Laboratories, Carlsbad, USA). For each sampling site one BFC and two NEC were included in the batch DNA extractions to monitor contamination and all eDNA extracts were stored at -20 °C for further analyses.

PCR amplification and the construction of the HTS libraries was done using the MiFish-U universal fish primers (Miya et al. 2015). These primers were selected based on the results of an in silico primer evaluation of the local fish diversity (Bylemans et al. 2018). Real-Time PCR reactions were run in triplicate for each sample with 25 µL individual reactions consisting of 0.20 μL of AmpliTaq Gold DNA Polymerase (5 U/μL; Applied Biosystems, Foster City, CA, USA), 2.50 μL of GeneAmp 10× Gold Buffer (Applied Biosystems), 2.00 μL of MgCl₂ (25 mmol/L; Applied Biosystems), 0.65 μL of GeneAmp dNTP Blend (10 mmol/L; Applied Biosystems), 0.20 μL UltraPure BSA (50 mg/ml; Invitrogen), 0.60 SYBR Green I Nucleic Acid Gel Stain (5X; Invitrogen), 1.00 μL of each primer (10 μmol/L), and 4.00 μL of template DNA and DEPC-treated water (Invitrogen). All PCRs were run on a Bio-Rad CFX96 Real-Time PCR System (Bio-Rad Laboratories, Hercules, USA) with PCR thermal cycling conditions consisting of 5 min at 95 °C; 45 3-step cycles of 30 s at 95 °C, 30 s at 61.5 °C, and 1 min at 72 °C; a final extension of 10 min at 72 °C and a melting curve with a stepwise increase of 0.1 °C/5 s from 60 to 95 °C. When negative control samples tested positive for fish eDNA, they were included in the HTS library preparation protocol. HTS libraries were constructed using a one-step Real-Time PCR amplification with fusion tagged primers (FTP) and the reaction and cycling conditions described above. Forward FTP consisted of the P5 sequencing adaptor, a custom forward sequencing primer, a 7 bp Multiplex Identification (MID-) tag and the MiFish-U forward primer. Reverse FTP contained the P7 sequencing adaptor, a 7 bp MID-tag and the MiFish-U reverse primer. MID-tags were generated using the edittag scripts and unique combinations of forward and reverse MID-tags were used to label amplicons from individual samples (Faircloth and Glenn 2012). For each sample, PCR reactions were performed in triplicate and amplicon libraries of 9 to 10 samples were pooled based on the average Ct-value. Pooled libraries were cleaned using Agencourt AMPure XP Beads (Beckman Coulter, Brea, USA) in a 1.2 volume ratio. Gel-electrophoresis (i.e. 2% agarose gel with a run time of 30 min at 120 V) was used to confirm the presence of a single amplicon and amplicons from each library pool were combined into a single super pool based on the band intensities. A total of 259 uniquely labelled amplicon libraries were combined in the final super pool that consisted of 60 libraries from the current study (i.e. 59 libraries from eDNA samples and 1 library from a BFC) and 199 libraries from other research projects. Unidirectional sequencing was conducted at the Ramaciotti Centre for Genomics (University of New South Wales) using the MiSeq v2 1x300bp sequencing kit.

Trimmomatic v.0.36 was used to trim technical sequences (i.e. sequencing adaptors and primers) from the sequencing reads (Bolger et al. 2014). Simultaneously low-quality bases (i.e. Q-score below 3) at the end of the sequencing reads were removed and a sliding window of 4-bases was used to trim reads when the average quality per base was below 15. Further bioinformatics filtering of the sequence reads followed the general workflow described in De Barba et al. (2014). The obitools scripts were used to assign sequence reads to their respective samples and remove sequences less than 150 bp in length and with an occurrence below 100 (Boyer et al. 2016). The threshold for removing low abundance sequences was determined experimentally so that all Actinopterygii sequences were removed from negative control samples. Sequences arising from PCR and sequencing errors were removed using the obiclean and obigrep scripts and taxonomic information was assigned to the sequence records using the ecotag script. The reference database used for taxonomic assignments was built using the standard vertebrate sequences from the EMBL data repository (release 132) and custom 12S sequences of all Actinopterygii species in the MDB (Bylemans et al. 2018).

Further filtering and analyses of the metabarcoding data was done using the packages tidyverse (Wickham 2016), inext (Hsieh et al. 2016) and vegan (Oksanen et al. 2007) in R version 3.4.1 (Suppl. material 2) (R Development Core Team 2010). The extra quality check revealed that some reads are likely to represent PCR and/or sequencing errors. Within the MR01 site Salmo trutta reads were highly abundant while the highest taxonomic assignment of a small proportion of the reads was the Salmo genus. The latter reads are likely to represent PCR and/or sequencing variants of S. trutta reads and were thus reassigned to S. trutta. A single sample of the MR05 site also returned reads assigned to Salmo salar and no other Salmonid species were detecting in this sample. As S. salar is unlikely to be present at this site these sequences likely represent contaminant DNA and were thus removed. As this site is a popular recreation spot one possible explanation could be the contamination by recreational users. Prior to further data analyses, species names were modified to ensure consistency between the expert opinion, electrofishing and metabarcoding data. As the barcoding region used in this study is unable to distinguish between species in the Galaxias genus all species were grouped into a single genus level taxonomy (Bylemans et al. 2018). Finally, all Hypseleotris species were renamed to H. klunzingeri as it is the only species of this genus to occur in the upper section of the Murrumbidgee River (Lintermans 2002).

The overall performance of the eDNA metabarcoding survey was evaluated by comparing the observed species richness with the estimated species richness obtained from both the expert opinion survey and the electrofishing data (Suppl. material 1).

To determine the impact of sampling intensity and spatial sampling replication the eDNA metabarcoding data was used to construct species accumulation curves (SAC), the total estimates species richness (S_est) and the number of samples required to detect 95% of S_est was determined. The data obtained from the MR01 site was excluded from the analyses as no variation in species detections were observed between samples. For the remaining sites, two sampling strategies were evaluated with the analyses using all available sample replicates or only those samples collected from the river-banks (i.e. LB and RB). The community data was first converted to presence/absence data and species accumulation curves (SAC) were constructed for each site and sampling strategy using the inext function (R package inext) (Hsieh et al. 2016). Additionally, the data was used to determine S_est (i.e. Chao2 estimates) and evaluate the number of samples required to detect 95% of S_est.

As sequencing depth (i.e. the number of reads per sample) could also influence the detection of species from metabarcoding analyses, the metabarcoding data with absolute read counts was used to determine the impact of sequencing depth. A custom R script relying on the R package vegan was used to first rarefy the community data to simulate different sequencing depths (i.e. 10,000; 20,000 and 40,000 reads per sample) (Oksanen et al. 2007). During this rarefaction both reads assigned to Actinopterygii species and reads discarded during the bioinformatics filtering process were considered (Suppl. material 1). Next, non-Actinopterygii reads were discarded and Actinopterygii reads with counts below 100 (i.e. the threshold used in the bioinformatics filtering process) were removed. The data was transformed to presence/absence data before using the inext function to construct SAC. Finally, rank abundance curves (RAC) were also constructed to express the relative read abundance for each species at each sampling site.

A permutational multivariate analysis of variance was used for a more in-depth evaluation of the impact of spatial sampling replication. Metabarcoding results were transformed to presence/absence data and community level differences between sampling transects and locations (i.e. left bank, mid-channel and right bank) were assessed using the adonis function within the R package vegan (Oksanen et al. 2007). Sampling transects and locations were set as independent variables and permutations were constrained within sites. When significant effects were observed, the simper function was used to evaluate the average contribution of each species to the overall dissimilarity.

After assigning the sequence reads to their respective samples an average of 31,365 sequence reads were obtained per sample with a minimum and maximum sequencing depth of 6,907 and 52,742 reads per sample, respectively. The overall quality of run was high (PhredQ30 score ≥ 91.17). The effect of the bioinformatics filtering processes on the number of sequencing reads for each sample is shown in the supplementary materials (Suppl. material 1).

The estimated/observed species richness for each sampling site and each survey method (i.e. expert opinion survey, electrofishing data and eDNA metabarcoding survey) is given in Figure 2. The results show the expected pattern of increasing richness with decreasing altitude. The total species richness of the two most upstream sites were highly similar for both the expert opinion and the eDNA metabarcoding surveys. However, the expert opinion data suggested that for each site one additional species was present but not detected during the metabarcoding survey (Galaxias sp. and Maccullochella macquariensis are believed to be present at the MR01 and MR02 sampling site respectively) (Suppl. material 1). The data from the three most downstream sites shows that the species richness obtained from the electrofishing data is consistently lower than the species richness derived from both the expert opinion and metabarcoding survey (Figure 2). Only on two occasions were species detected by the electrofishing surveys and not by the eDNA metabarcoding survey (i.e. Macquaria australasica at the MR03 site and Macquaria ambigua at the MR04 site) (Suppl. material 1). Contrastingly, three, five, and six species were only detected by the eDNA metabarcoding survey at the MR03, MR04 and MR05 site respectively (Suppl. material 1).

Figure 2.

The species richness for each sampling site within the Murrumbidgee River obtained from the standard electrofishing surveys, the expert opinion survey and the eDNA metabarcoding survey.

The SAC constructed for each site using all samples and the river-bank samples only show that the number of samples needed to characterize the species community generally increases with decreasing altitude and an increase in species richness (Figure 3A). For both the MR02 and MR03 sites sampling the river-banks only is an appropriate strategy with 12 and 13 samples being sufficient to detect ≥ 95% of S_est, respectively (Figure 3A). The results of the MR04 site reveal that approximately 16 samples collected from both river-banks would be adequate to detect 95% of S_est. However, increasing the spatial sampling replication reduced the number of samples required to characterize the species richness (i.e. the collection of approximately 13 samples from both river-banks and the mid-channel allows for the detection of 95% of S_est) (Figure 3A). Finally, at the most downstream site (i.e. MR05) the results show that only sampling the river-banks will result in an underestimation of the true biodiversity and collecting approximately 15 samples from both river-banks and mid-channel is sufficient to detect ≥ 95% of S_est (Figure 3A).

Figure 3.

Species accumulation curves (SAC) (A–B) and rank abundance curves (RAC) (C) for the four most downstream sampling sites in the Murrumbidgee River. SAC were constructed to compare two different sampling strategies (A) (i.e. using all available samples and only samples collected from the river-banks) and different sequencing depths (B). The vertical dashed lines in panel A show the number of samples needed to detect ≥ 95% of the estimated species richness (i.e. based on the Chao2 estimates) for each site and sampling strategy.

The results from the SAC derived from the data with different simulated sequencing depths show that, for most sites, increasing sequencing depth only moderately improves species detections (Figure 3B). While there is no obvious effect of species diversity on the impact of sequencing depth, when considering both the SAC and the RAC the results show that sequencing depth only has an impact for the MR02 and MR05 sites (Figure 3B, C). The results of the RAC show that for these sites some species have a very low proportional abundance of sequence reads (i.e. below 1% of the total sequence reads per samples) and these species are likely to be missed when using a sequencing depth of 10,000 reads per sample.

The results of the analysis of variance of the community data reveal the need for spatial sampling replication at low altitude sites with a high species richness and is consistent with the results of the SAC. The analysis of variance of the community data showed that sampling location had a significant effect on the community data (R² = 0.01931, p-value = 0.043) while sampling transect did not (R² = 0.00618, p-value = 0.876). The results show that the overall community dissimilarity between sampling locations is relatively low for the MR02 site compared to the three most downstream sites (Figure 4A). No clear patterns were observed in the overall community dissimilarity for the MR03 and MR04 sampling sites. The results from the MR05 sampling site, however, show that the community dissimilarity is higher between the mid-channel sampling location and the river-bank sampling locations (i.e. LB and RB) (Figure 4A). In all three downstream sampling locations the detections of Galaxias sp. appear to be relatively important in determining the community dissimilarity between mid-channel samples and river-bank samples (Figure 4B). Galaxias sp. detections were overall higher for the river-bank samples collected at the MR03 and MR05 sites, while the reverse pattern was found for the MR04 site (Suppl. material 1). The most noticeable community differences between sampling locations were found in the MR05 sampling site. Here four fish species (i.e. Gambusia holbrooki, Galaxias sp., Maccullochella peelii and S. trutta) contribute to the community dissimilarity between the mid-channel sampling location and the river-bank sampling locations but did not contribute to the dissimilarity between the two river-bank locations (Figure 4B). While G. holbrooki and S. trutta were exclusively detected in the mid-channel samples collected at the MR05 site; M. macquariensis and Perca fluviatilis detections were only found in the river-bank samples (Suppl. material 1).

Figure 4.

Overall dissimilarity between the fish community data obtained from the samples collected from the left bank (LB), mid-channel (MC) and right bank (RB) for the different sampling sites (A). The heat map shows the average contribution of each species to the overall community dissimilarity (B).