eDNA along Houdong riverine zonation in Taiwan

Study Extent

This is a one-time sampling event. Water samples were collected along the Houdong River in Jiaoxi Township, Yilan County, Taiwan. The Houdong River is a popular tourist destination running across the city of Yilan. The river system originates east of the Sidu mountains and flows through primary and secondary forests, agricultural lands (rice), and developed areas (light industrial and residential) until it eventually drains into the Pacific Ocean. We selected four sites along the Houdong River. We collected water samples at each site for eDNA analysis and measured in-situ water quality.

Sampling

Water samples for eDNA were collected from the near surface of the river. Before the collection, the water containers were rinsed with the local water at each sampling site. 1L of water was collected for eDNA analysis, and 200mL of water was collected for water quality measurements and were measured using multiple handheld probes on-site at each sampling site. Temperature, pH, and dissolved oxygen were measured through a multi-parameter meter (Multiline® Multi 3620 IDS, WTW, Weilheim, Germany) equipped with an IDS pH electrode (SenTix 940, WTW) and an optical IDS dissolved oxygen sensor (FDO® 925, WTW). Turbidity was measured by the turbidity meter (TUB-430, EZDO, Taiwan). Salinity was determined by the sanity refractometer (2491 MASTER-S/Milla Salinity Refractometer, Atago, Japan). Three replicate measurements of each parameter were taken. Water samples were transported and processed through two filtration steps at the Marine Research Station (MRS, Yilan, Taiwan) of the Institute of Cellular and Organismic Biology (Academia Sinica, Taipei, Taiwan). In the first step, the larger particles were removed by a 75µm pore size sieve. Afterward, the 1L water sample from each site was vacuum-filtered with a 0.22µm filter (PC651-0024, GeneDireX, USA). The filter membranes were placed in sterile Petri dishes and stored at -80°C until DNA extraction.

Method steps

1. Wet lab process: The DNA was extracted at the Biodiversity Research Center (Academia Sinica, Taipei, Taiwan). Each filtered membrane was cut into quarters. 3 of 4 pieces of filtered membranes were used in the study as 3 experimental replicates. The final quarter was saved as the sample backup. DNA from each quarter membrane piece was extracted using the Presto™ Stool DNA Extraction Kit (STLD100, Geneaid Biotech Ltd., Taiwan) following the manufacturer's instructions (Instruction Manual Ver. 10.21.17). The quality and quantity of the extracted DNA was assessed using a Nanodrop 2000 (Thermo Fisher Scientific Inc., USA) and the Qubit 4 dsDNA High Sensitivity Assay Kit (Thermo Fisher Scientific Inc., USA). The MinibarF1 (5'TCCACTAATCACAARGATATTGGTAC) and MinibarR1 (5'GAAAATCATAATGAAGGCATGAGC) primers that were designed by Meusnier et al. (2008), were used to amplify the 5' region (ca. 120-150bp) of the mitochondrial Cytochrome c oxidase I (COI) gene. The universality of the primers was recommended for distinguishing the highly diverse DNA from the environmental mixture. We conducted PCR using a one-step single-indexed approach, with a 13bp tag attached to the MinibarR1 primer. The PCR reaction volume was 16μl, which included 8μl KAPA HiFi HotStart ReadyMix (KK2602, Roche Molecular Systems Inc., USA), 5μl ddH20, 1μl of each primer (10μM), and 1μl of DNA template. To optimize the protocol, we performed a preliminary PCR using an annealing temperature gradient and found that 54°C gave the best results. The PCR mixture was denatured at 95°C for 15 minutes, followed by 35 cycles of 94°C for 30 seconds, 54°C for 30 seconds, and a final elongation at 72°C for 10 minutes. The PCR products were checked on a 1.5% agarose gel and quantified with the Invitrogen Qubit 4 fluorometer (Thermo Fisher Scientific Inc., USA). Afterwards, all the PCR products were pooled in one tube for next-generation sequencing. Sequencing was performed on the Illumina NovaSeq 6000 platform with 2*150 paired-end reads by Genomics Co., Taipei, Taiwan.
2. Data processing and analyzing: The Illumina raw reads were demultiplexed by Genomics Co., Taipei, Taiwan. FastQC (v0.11.9; https://github.com/s-andrews/FastQC) was used to the quality-checked. The forward and reverse primers of the demultiplexed reads were trimmed using Cutadapt (version 4.2; Martin 2011). The USEARCH platform (v11.0.667; Edgar 2010) was used to verify if the primer sequences were completely removed from the demultiplexed reads. The "denoised-paired" function was used to create an amplicon sequence variant (ASV) data output from the demultiplexed reads. The "denoise-paired" function in QIIME2 (v2023.2.0; Bolyen et al. 2019) can automatically trim, filter, denoise, merge reads, and remove chimeric reads in one step. The maximum expected error of forward and reverse reads was set to 1.0. The reads with a quality score of less than 20 were truncated. The minimum overlap length for the forward and reverse reads merger was set to 16 bp. Other parameters followed the default settings in the "denoise-paired" function. No reads were trimmed or truncated during the ASV creation process. The ASV output was then further clustered through the "cluster-features-de-novo" function provided by QIIME 2. ASVs with more than 97% sequence identity were clustered into one operational taxonomic unit (OTU). The sequences that were shorter than 100bp in the ASV and OTU results were discarded. Taxonomic assignments were conducted using Constax (v2.0.18; Liber et al. 2021) against the MIDORI database (vGB250; Machida et al. 2017). The R package phyloseq was used to analyze the preprocessed sequencing data (v1.40.0; McMurdie and Holmes 2013). The pie chart figure was produced with ggplot2 (v3.4.2; Villanueva and Chen 2019). Lowest taxon level-annotation of each OTU was extracted to perform a secondary species mapping to the GBIF Backbone Taxonomy (GBIF Secretariat 2011) using the "name_backbone_checklist" function in R package rgbif (v3.7.7; Chamberlain et al. 2022) in R (R Core Team 2023).
3. Open data and code: DNA sequence data have been deposited on ENA at EMBL-EBI under accession number PRJEB60905. We converted the occurrence data into Darwin Core Archive standard (Darwin Core Task Group 2009) and validated data in the event and occurrence data sheets using the GBIF Data Validator (Global Biodiversity Information Facility 2017). We then published the dataset containing one event core and four extensions following GBIF registered extensions: occurrence, DNA-derived data, measurement or facts, and resource relationship using the Integrated Publishing Toolkit (IPT) of GBIF installed under the Taiwan Biodiversity Information Facility (TaiBIF). All source code used in the project can be found in the project's GitHub repository.

We detected eukaryotic organisms in the water samples using the COI mitochondrial gene. A total of 2,736 OTUs were identified and 421 of the OTUs were assigned to at least the kingdom level.

1. Eukaryota
   rank: domain
2. Animalia
   rank: kingdom
3. Chromista
   rank: kingdom
4. Fungi
   rank: kingdom
5. Plantae
   rank: kingdom
6. Protozoa
   rank: kingdom
7. Amoebozoa
   rank: phylum
8. Bryophyta
   rank: phylum
9. Cryptophyta
   rank: phylum
10. Nemertea
    rank: phylum
11. Sulcozoa
    rank: phylum
12. Annelida
    rank: phylum
13. Bryozoa
    rank: phylum
14. Gastrotricha
    rank: phylum
15. Ochrophyta
    rank: phylum
16. Tracheophyta
    rank: phylum
17. Arthropoda
    rank: phylum
18. Chaetognatha
    rank: phylum
19. Glomeromycota
    rank: phylum
20. Oomycota
    rank: phylum
21. Zygomycota
    rank: phylum
22. Ascomycota
    rank: phylum
23. Chlorophyta
    rank: phylum
24. Haptophyta
    rank: phylum
25. Platyhelminthes
    rank: phylum
26. Basidiomycota
    rank: phylum
27. Chordata
    rank: phylum
28. Mollusca
    rank: phylum
29. Rhodophyta
    rank: phylum
30. Blastocladiomycota
    rank: phylum
31. Cnidaria
    rank: phylum
32. Mycetozoa
    rank: phylum
33. Rotifera
    rank: phylum

Houdong River (猴洞坑), Jiaoxi Township, Yilan County, Taiwan