Assessing freshwater fish biodiversity of Kumbe River, Papua (Indonesia) through environmental DNA metabarcoding

Keywords: biodiversity, eDNA, freshwater fish, Indonesia, metabarcoding, Papua, river, tropical.

Defined by Walter Rosen in 1986 during the first American forum on biological diversity, biodiversity represents the variety and variability among living organisms and ecosystems where they occur (Pereira et al. 2013; Bartkowski et al. 2015). The challenge for biodiversity conservation is to understand when and how major perturbations such as climate change cause significant impact (Yachi and Loreau 1999; Mace et al. 2012). Thus, it is the responsibility of individual countries to preserve biodiversity for its environmental, productive, consumptive, social, ethical and aesthetic value (Pauchard 2017).

Indonesia is one of the most biodiverse countries, and its freshwater fish biodiversity is no exception. A total of 1200 native fish species have been reported across the archipelago (Kottelat and Whitten 1996; Hubert et al. 2015). The inventory of Indonesian freshwater species is ongoing, and recent studies have suggested that fish diversity may be underestimated (Kadarusman et al. 2012; Hubert et al. 2019; Sholihah et al. 2020, 2021a, 2021b). Thus, the incomplete understanding of fish biodiversity presents challenges for effective conservation and management. While traditional methods can be implemented to document biodiversity, these can be cost-prohibitive and time consuming. In addition, it also requires detailed expertise in fish taxonomy, which becomes limiting as the number of taxonomists declines. Optimising monitoring strategies will be critical to determine change in ecosystems; and monitoring over time is a crucial step toward sustainable use of resources (Xiao et al. 2016; Trebitz et al. 2017). This challenge is particularly pronounced in Indonesia, where aquatic biodiversity-related research output is not balanced with biodiversity complexity (Wibowo et al. 2018; Gustiano et al. 2021; Kurniawan et al. 2021) and understanding of the dynamics behind the origin and maintenance of freshwater biodiversity is still fragmentary (de Bruyn et al. 2013, 2014; Sholihah et al. 2021a, 2021b).

Threats to Indonesian aquatic biodiversity have escalated in recent decades through agricultural development (Cleary and DeVantier 2011; Imai et al. 2018; Austin et al. 2019), introduction of exotic species (Herder et al. 2012; Dahruddin et al. 2017), and pollution (Garg et al. 2018). As such, locally decreasing ichthyodiversity is of concern, and taxonomic knowledge gaps are challenging (Dahruddin et al. 2017; Hubert et al. 2019). Papua, the western most Indonesian province of the Island of New Guinea, contains 16% of known Indonesian freshwater fish diversity and 20% of endemic species, including all fish species within the genus Melanotaenia (Rainbowfish) (Hubert et al. 2015). New species or expanded distributions of known species are still being reported (Nugraha et al. 2015; Wibowo et al. 2017; Ditya et al. 2018). However, much of the fish biodiversity of Papua is still undiscovered (Kadarusman et al. 2012; Koh et al. 2013; Hubert et al. 2015).

An alternative approach to identify fish biodiversity is to collect DNA shed by the target species into the water. All aquatic organisms release DNA into the environment via mucus, urine, faeces, or dead tissue (Goldberg et al. 2016; Carraro et al. 2018), where it is then known as environmental DNA (eDNA). Water samples containing eDNA from the entire aquatic community (e.g. bacteria, algae, fish) can be filtered to extract eDNA and the target species can be identified using the relevant primer sets (Turner et al. 2015; Pont et al. 2018). Environmental DNA sampling can be used to detect a single species or multiple species (known as metabarcoding) (Taberlet et al. 2012; Rees et al. 2014; Pont et al. 2018). Although species-specific eDNA sampling may be more sensitive than metabarcoding – especially for rare species (Bylemans et al. 2019) – eDNA metabarcoding is a rapid, cost-effective, and non-invasive technique (Janosik and Johnston 2015; Seymour 2019; Sigsgaard et al. 2020). Environmental DNA has been successfully implemented in the surveillance of rare, invasive, or migratory fish species in various aquatic systems (Nathan et al. 2015; Sigsgaard et al. 2015; Simpfendorfer et al. 2016; Rice et al. 2018; Itakura et al. 2019; McElroy et al. 2020). However, eDNA metabarcoding has rarely been applied in Indonesia, with only three studies in freshwater environments. For instance, eDNA was successfully implemented in detection of invasive alligator gar Atractosteus spatula (Ulayya et al. 2020) and crayfish Cherax quadricarinatus (Djalil et al. 2018) in lakes of West Java. Assessing fish biodiversity using eDNA was carried out in the Maninjau lake of Sumatra, where 26 species were detected in the lake, including five native and 21 exotic fish species (Roesma et al. 2021).

The objective of this study was to evaluate the utility of eDNA metabarcoding to rapidly characterise the fish community of the Kumbe River in West Papua and assess its potential in a tropical context as a suitable alternative method to traditional sampling to characterise fish diversity.

Sampling was undertaken at four locations in the Kumbe River, which is located in the southern part of Papua Island (Fig. 1). The Kumbe River is a meandering tropical river with a total length of 242 km and width ranging from 97 to 700 m (Anonymous 2010). The average annual rainfall of this area is 2160 mm and average annual temperature of 23°C (Peel et al. 2007). We collected seven 15 mL water samples (a total of 28 samples) from each of four sampling locations (Fig. 1). Immediately after collection, we added 1.5 mL of 3 M sodium acetate and 33 mL of absolute ethanol, mixed with the water sample and then stored it at room temperature in the field. Once in the laboratory, water samples were placed at −20°C for DNA preservation until DNA extraction.

Fig. 1.  (a) Map of the study location in the Kumbe River, Papua Island, Indonesia and photos of each sampling site, (b) Alfasera: S′07°24′23.9″; E′140°37′04.6″, (c) Inggun: S′07°59′05.3″; E′140°27′53.3″, (d) Yakui: S′08°02′26.3″; E′140°31′47.5″, and (e) Sakor: S′08°06′51.6″; E′140°29′98.9″.

Water samples were centrifuged for 35 min at 5500g at 24°C to precipitate DNA and suspended material, and supernatant was discarded (Valiere and Taberlet 2000). Precipitates were resuspended in 100 μL ultrapure water and DNA was extracted using the QIAamp Tissue Extraction Kit (Qiagen). We extracted DNA from each precipitate individually and DNA extraction was performed in a DNA laminar flow hood to avoid contamination. A negative control containing all reagents but no eDNA was included in each batch of extractions to monitor contamination. The concentration of extracted DNA ranged from 1 to 20 ng/μL. Polymerase chain reaction (PCR) amplification was performed using the primers 5′-TGCCAACGGAGCATCATTC-3′ and 5′-ATAAAGGTAGGAGCCGTAGT-3′, which amplify a 79 base pairs (bp) segment of the mitochondrial Cytochrome b (Cyt-b) (Ficetola et al. 2008). Each PCR reaction contained 2.5 μL of 1 × KAPA HiFi Hotstart ReadyMix (Kapa Biosystems), 0.5 μL of each primer (1 mM) and 5 μL of DNA extract; 4.0 μL PCR grade water was added up to a reaction volume of 12.5 μL. A negative control was included with each PCR batch to monitor the contamination of the reagents used. The PCR conditions were as follows: one cycle at 95°C for 5 min; 40 cycles at 95°C for 30 s, 50°C for 30 s, 72°C for 30 s; and one cycle at 72°C for 5 min. PCR products were visualised under UV light using a 2% agarose gel stained with Midori Green Advanced (Nippon Genetics Europe). PCR products were purified using the NucleoMag 96 PCR Kit (Macherey-Nagel).

We pooled individually tagged DNA samples into two different sequencing libraries. Libraries were set up as described by Meyer and Kircher (2010) but modified as follows for Ion Torrent Personal Genome Machine (PGM) sequencing. The adapter mix was prepared using an adapter marked with a multiplex rating (MID). Adapters are designed to include individual MIDs and to match Ion Torrent-specific priming sites. The total molar concentration of the adapter mixture was 20 M for each adapter. We used a 100 bp long PCR product as the positive control DNA template. We carried out blunt-end reaction with half of the volume, i.e. we added 10 μL of blunt-end master mix to 25 μL of the sample.

The reaction was purified using solid phase reversible immobilisation (SPRI) beads with a template: bead ratio of 1:1.8. In the post-ligation purification steps, the template:bead ratio was 1:1. In order to verify adapter ligation success, the SPRI-purified MID-tagged positive control library and a subset of MID-tagged library aliquots were separated side by side with non-MID-tagged templates on 2% agarose gel stained with Midori Green Advanced electrophoresis (45 min, 95 V) and visualised using UV light. We measured DNA concentration of the samples using a Qubit Fluorometer (Life Technologies) and mixed amplicon libraries subsequently in equimolar ratios. The stock library was re-amplified in two different reactions as follows: 5 μL pooled library stock was added to a master mix comprising of 5 U of Herculase II polymerase (Agilent technologies, catalogue number 600677), 1 × Herculase II reaction buffer, 25 mM each deoxynucleoside triphosphate, 10 μM each primer, and added PCR grade water up to 50 μL. We designed the re-amplification primers based on the sequence of the IonTorrent-specific priming sites and thus generated millions of MID-tagged copies including binding sites necessary for subsequent sequencing with IonTorrent technology. The thermocycling profile included a 30 s denaturation at 98°C followed by 15 cycles consisting of a 20 s denaturation at 98°C, a 30 s annealing at 64°C, and a 30 s elongation at 72°C. Final elongation was conducted at 72°C for 5 min. To clean the re-amplified library, size-selection was conducted by separating the entire library using 2% Size-Select Agarose E- Gel and E-Gel Electrophoresis System (Life Technologies) according to the manufacturer’s instructions. The DNA concentration of amplified library pool was measured using a Qubit Fluorometer. The library pool stock was then diluted to a final concentration of 26 pmol. For template preparation, an 18 μL aliquot of the library dilution (approximately 2.8 × 10⁸ molecules) was transferred into the sequencing reaction set-up. Emulsion PCR and Ion Torrent PGM sequencing were carried out on two 314 chips (Life Sciences, catalogue number 4462923) according to the manufacture’s protocol (Publication Part Number: 4471974 Rev. C).

The resulting reads were binned and renamed by MID, i.e. the original individual samples using the software Geneious Pro 6.1. Then the original-specific primers (Ficetola et al. 2008) were trimmed off, reads were trimmed for poor quality parts using a 0.05 error probability limit and reads shorter than 100 bp were excluded using Geneious Pro 6.1. Thus, only reads comprising both the Ion Torrent adapter with a MID and the original primer were passed on for further analysis. Subsequent analyses were performed employing supercomputers at the IT Centre for fisheries (Yokohama, Japan). Sequences were assigned to species using the BLASTN 2.2.25+ algorithm against the reference database (GenBank). We examined Basic Local Alignment Search Tool (BLAST) results showing a sequence similarity of at least 94% and we also considered species distribution when confirming species identity (Wibowo et al. 2017).

The BLAST output was imported into the program MEGAN (Huson et al. 2007) to further summarise the results. We employed default parameters, except for following settings: Min support = 1, Min score = 100, Min complexity = 0. We utilised a built-in comparison tool in MEGAN for database comparison. We discarded all hits below 94% identity. We also discarded hits with e-value over 1e-20. We further analysed the hits one by one to see if there were reads matching more than one species in GenBank. In case of multiple hits, we (a) identified the hit with the longest alignment length, (b) identified the hit with a match to an adult specimen over a juvenile specimen (given the greater confidence in morphological identification of an adult specimen), and (c) reviewed the database entry for errors and omitted incorrect entries.

A total of 629 976 raw reads were detected by sequencing of amplicons from 28 water samples and filtered to 297 438 reads. Thus, altogether 332 538 reads (52.79%) were discarded from further analysis after quality control. The numbers of filtered sequences obtained from Alfasera, Yakui, Sakor, and Inggun were 71 627, 89 014, 81 030 and 55 767 reads, respectively. BLAST queries resulted in sequence similarities between reads and reference sequences in Genbank ranging between 94 and 98% (Table 1). A total of 23 freshwater fish species were detected (Table 1).

Table 1.  The list of fish species detected at Kumbe River collected based on eDNA and previous conventional sampling.

Species diversity was greatest at Inggun (21 species) and Yakui (21 species), followed by Sakor (19 species) and Alfasera (16 species) (Fig. 2a–d). The number of reads per species varied considerably among the seven sampling replicates within sites, with reads proportion ranging from 0.07 to 0.27 for Oreochromis niloticus, at Alfasera for instance (Fig. 2a). Some rarer species were not represented in all sampling replicates within sites including Trichopodus microlepis, Arius arius and Smilosicyopus bitaeniatus, which were not represented in replicates one and five at Yakui, and Melanotaenia goldiei, which was not represented in replicate six at Inggun.

Fig. 2.  Relative abundance of species-specific reads across PCR replicates at (a) Alfasera, (b), Yakui, (c) Sakor, and (d) Inggun.

Fourteen of the 23 species were detected at all sites (M. goldiei, O. niloticus, Craterocephalus randi, Craterocephalus stercusmuscarum, Craterocephalus nouhuysi, Craterocephalus gloveri, S. bitaeniatus, Neogobius melanostomus, Oryzias haugiangensis, Atherinella schultzi, Acanthopagrus berda, Diplospinus multistriatus, Anguilla australis and Labeobarbus natalensis), five species were detected at three sites (Melanotaenia albimarginata, M. sp, Neoarius graeffei, Neoarius berneyi and Oryzias latipes), two species were detected at in two sites (Oxyeleotris lineolata and Pseudanthias dispar) and two species were detected at one site (A. arius and T. microlepis) (Table 1, Fig. 3). A total of 16 species were native (M. goldiei, M. albimarginata, M. sp., C. randi, C. stercusmuscarum, C. nouhuysi, C. gloveri, N. melanostomus, N. berneyi, N. graeffei, A. berda, A. australis, D. multistriatus, O. lineolata, P. dispar and A. arius) and seven species were exotic (O. niloticus, S. bitaeniatus, O. haugiangensis, O. latipes, A. schultzi, L. natalensis and T. microlepis) (Table 1). If reads were ranked by abundance (Table 1), three species dominated all sites, consisting of M. goldiei (32.03% of the all reads), C. randi (31.19% of all reads), and O. niloticus (28.77% of reads). The 23 detected species belonged to 13 families and 16 genera. Atherinidae (35% of all sites; four species) was the most abundant family in terms of sequencing reads, followed by Melanotaeniidae (32% of all sites; three species) and Cichlidae (29% of all sites; two species). Sequences from these three families dominated all studied locations (Fig. 4a). These families mostly represented the genera Craterocephalus, Melanotaenia and Oreochromis (Fig. 4b).

Fig. 3.  Venn diagram of identified number of shared fish species among four sampling locations.

Fig. 4.  Fish community structure (proportion of total reads) at family (a) and genus (b) levels detected using eDNA metabarcoding at four different sites of the Kumbe River.

A total of 23 species have previously been detected in the Kumbe River using conventional methods and DNA barcoding (Table 1). We detected five of these species (N. graeffei, M. goldie, C. stercusmuscarum, O. niloticus and A. arius) using metabarcoding, and an additional 18 species (69.56% native species) not previously recorded from the Kumbe River.

Environmental DNA was successfully deployed to characterise freshwater fishes species richness and species occurrences at the Kumbe River. This represents one of the first eDNA metabarcoding assessments of fish biodiversity in Indonesia. Interestingly, N. graeffei or blue catfish was the dominant fish species (90%) in the Kumbe River based on monthly fishing over 2 years (Ditya et al. 2018). However, eDNA metabarcoding revealed only five N. graeffei reads from the whole dataset, while 91.9% of reads were mapped to only three species (M. goldiei, C. randi and O. niloticus). Thus, the lack of consensus between earlier studies employing traditional fishing methods (Wibowo et al. 2017; Ditya et al. 2018) and current eDNA analysis requires further investigation.

Our results showed that only five of the species we detected had been previously recorded (Wibowo et al. 2017; Ditya et al. 2018), and 18 represented new records, highlighting the limited knowledge of fish biodiversity in the Kumbe River. The detection of 18 new species is an important addition and highlights the critical role metabarcoding could have in rapidly describing fish biodiversity in tropical ecosystems. However, given that very low number of mapped reads (<50) for most of the species, the occurrence of these fishes in the River Kumbe needs to be validated by traditional approaches. Furthermore, it is possible that fish biodiversity in the Kumbe River is greater than reported here given many sequences were produced in this study that could not be identified due to gaps in the reference database. This demonstrates that it is critical that the reference databases be expanded to include more taxonomic groups to improve our understanding of tropical fish biodiversity (Schenekar et al. 2020) by taking every opportunity to collect fin clips from species where identification has been verified.

The inability of metabarcoding to detect 17 species previously recorded in the Kumbe River could be due to a range of factors. From a technical perspective, several reasons may account for this discrepancy: (1) some species may not have been present at the time of eDNA sampling (e.g. migratory species); (2) failure to detect a species using eDNA may be due to low abundance, for example, earlier studies have shown that eDNA metabarcoding may generate false negative results for rare species or lineages (Sato et al. 2017; Burian et al. 2021; Xiong et al. 2022); (3) the incomplete detection of species based on eDNA may be caused by insufficient volume of water filtered, filtration method, PCR inhibitors, primer bias, preferential amplification of abundant species, or an incomplete database of reference sequences (Deiner et al. 2017; Majaneva et al. 2018; Bessey et al. 2020; Schenekar et al. 2020; Rojahn et al. 2021).

We chose to use the Cyt-b gene fragment as it is commonly applied for fish metabarcoding studies. Some studies have demonstrated that Cyt-b may not perform as well as other regions, such as the 12S rRNA gene (Zhang et al. 2020; Shu et al. 2021). However, not all studies show consistent patterns and there may be distinct advantages and disadvantages with both markers (Hänfling et al. 2016). It is possible that the number of species detected may be maximised using multiple markers, which also enables the verification of taxa detected in each sample (Lecaudey et al. 2019). We therefore recommend that future metabarcoding studies in Indonesia evaluate the performance of multiple mitochondrial and nuclear markers to identify the most powerful targets for species identification.

From an ecological perspective, there are several limitations to characterising biodiversity based on eDNA metabarcoding. Environmental DNA production, degradation and transport are the three key factors that determine the quantity of eDNA in the water (Hansen et al. 2018). Thus several abiotic factors such as temperature, pH, dissolved oxygen and organic matter can degrade eDNA and reduce species detectability. In this study we did not account for any environmental factors, and thus detectability of some species may have been reduced. If eDNA metabarcoding was to become more widely applied in the tropical waters of Indonesia for conservation purposes – particularly for detection and quantification of specific species – we recommend further research into the role that abiotic and biotic factors play in eDNA production and degradation (Stewart 2019). In addition, false-positives detections may be caused by misidentification of taxa or contamination during field sampling or in the laboratory (Sato et al. 2017; Fujii et al. 2019). Thus, more comparative studies at species-rich regions are needed to evaluate when eDNA metabarcoding estimates of species richness are on par with, or potentially better than traditional fisheries methods (Ficetola et al. 2015; Olds et al. 2016; Cilleros et al. 2019). Nevertheless, despite the limitations, we have demonstrated that eDNA metabarcoding is a promising method to assess biodiversity in Indonesia. We recommend that eDNA metabarcoding be deployed concurrently with traditional methods for optimal results while methods are being optimised and reference databases are being improved.

Fish biomonitoring using eDNA has been implemented in the detection of invasive species and assessing their distribution (Nathan et al. 2015; Dunker et al. 2016). The invasive Nile tilapia O. niloticus was detected at all sites, which is concerning given it is implicated in extirpations of native fish species (Ogutu-Ohwayo 1990; Ligtvoet et al. 1991). The Nile tilapia is an important source of animal protein across its native range, and as such, it has been widely introduced in Papua for aquaculture purposes. Environmental DNA metabarcoding in Indonesian waterways may be a rapid and cost-effective method to track new introductions or expansion of Nile tilapia.

The Kumbe River drains into the Arafura Sea and is expected to contain a diverse fish community comprised of potamodromous species as well as species that require access to both freshwater and seawater to complete their lifecycle. Given the Kumbe River is a pristine aquatic habitat where land use has had minimal negative impacts on the environment (Wibowo et al. 2017; Lasmana et al. 2018), eDNA metabarcoding could be used to determine if migratory species are able to access upstream habitat as expected (Duda et al. 2021). In addition, eDNA metabarcoding could also be used to monitor spawning patterns of migratory fish species (either potamodromous or diadromous) (Yamanaka and Minamoto 2016; Thalinger et al. 2019). For example, if migratory species diversity is similar at upstream and downstream locations, this indicates good habitat connectivity between studied locations (Yamanaka and Minamoto 2016; Carraro et al. 2018). We detected A. australis, known as catadromous species, at the most upstream site (Alfasera). This likely suggests that the lower reaches of the Kumbe River, where this study was conducted, have sufficient connectivity to allow fish passage upstream. However, as this species is able to climb over migration barriers (Jellyman et al. 2017), A. australis presence upstream may not indicate that there is clear passage upstream for other species with lower upstream migration capabilities.

Environmental DNA metabarcoding shows great promise as a novel tool for fish biomonitoring programs in tropical freshwater habitats. It is sensitive, non-invasive and cost-effective as well as enabling a rapid assessment for detection of invasive, rare or endangered species. It provides an alternative approach when traditional sampling is time consuming and requires taxonomic expertise. However, our results highlight that eDNA metabarcoding cannot be routinely applied without further research into sampling design, marker choice, seasonal data and an improved reference databases to reduce false negative detections. Despite this, eDNA metabarcoding will play an important role in the future characterisation of the fish community in tropical river systems.

The data used in this manuscript are organised by the first author (Arif Wibowo). Access to this data can be negotiated with him.

The authors declare no conflicts of interest.

This study was supported financially by the Research Institute for Inland Fisheries and Extensions (RIIFE), Division of Genetics, Department of Biology, Turku University, the Erasmus Mundus Research Fellowship Programme, Academy of Finland, and the Estonian Ministry of Education and Research (Institutional Research Funding Project IUT8–2).