Detection of terrestrial mammals using environmental DNA during heavy rainfall events and associated influencing factors

Study area and target species

Portions of this text were previously published as part of a preprint (https://doi.org/10.22541/au.174312906.60655689/v1). The study area included three rivers within the Kiyotake River catchment in Miyazaki Prefecture, southwestern Japan: the Mizunashi, Oka, and Tagami Rivers as well as the irrigation canal connected to the latter (Fig. 1A). The Kiyotake River is 28.8 km long and its area covers 166.4 km² (data: https://www.pref.miyazaki.lg.jp). As land use influences the intensity of surface runoff (Guzha et al., 2018), we utilized the Advanced Land Observing Satellite (10-m accuracy, data available: https://earth.jaxa.jp/ja/data) to classify land use into 12 categories (e.g., Urban, Field, Deciduous Broadleaf Trees, Deciduous Needleleaf Trees, Fig. 1B), and further classified into four categories: forest (e.g., Deciduous Broadleaf Trees, Deciduous Needleleaf Trees), farmland (Paddy Field, and Field), urban areas (Urban), and others (e.g., Water Area, Grassland). The proportion of each land use was derived by dividing the number of meshes corresponding to the land-use type (e.g., forested area) by the total number of meshes in the watershed of the studied sites. The major watershed land used in the studied rivers are forest in the Mizunashi River (71.57%), forest and farmland in the Oka River (39.25% and 31.25%, respectively), and farmland and urban areas in the Tagami River (43.24% and 21.38%, respectively).

Our study focused on livestock cattle species, B. taurus. As a region with a well-developed livestock industry, Miyazaki City contains multiple cattle sheds, providing widespread sources of eDNA release. Due to the small scale of each cattle shed (200∼300 m²), we hypothesize that the abundance of this model species at each site was comparatively low and that it is confined within enclosures and designated areas. Each cattle shed is situated at a certain distance and elevation from the corresponding river, making it physically impossible for the cattle to enter these natural water bodies. Therefore, the direct use of rivers by target species is negligible.

Local regulations prohibit wastewater from being directly discharged into rivers, thereby limiting eDNA from wastewater sources. The sewer system is a combined system that treats wastewater at water treatment facilities before it is discharged. This species has a low risk of false positive results due to the absence of individuals from the same species in the wild. At least one cattle shed raising individuals outdoors is present in each river catchment (Fig. 1). The specific geographic coordinates of the sampling locations were: the Mizunashi River (31.829925°N, 131.362669°E), the Oka River (31.85162°N, 131.386959°E), the Tagami River (31.838143°N, 131.417732°E), and the irrigation canal (31.839362°N, 131.414762°E) and geographic coordinates of the cattle sheds were: the Mizunashi River (31.828235° N, 131.359762°E), the Oka River (31.852260°N, 131.385643°E) the Tagami River and the irrigation canal (31.839116°N, 131.413398°E).

Sampling

A scheme of field investigations and experiments is presented in Fig. 2. Prior to water sampling, one and 0.25 L plastic bottles, tethered buckets, cooler boxes, and a scoop were sterilized with 10% bleach for at least 30 min to prevent contamination by DNA and microorganisms other than the target species. We collected samples multiple times after and during moderate to heavy rainfall events from 2021 to 2023. Specifically, we collected samples twice a day (Sampling 1 and 2) on September 14, 2021, November 22, 2021, May 12, 2022, and June 2, 2023, to understand the temporal pattern of eDNA detection. The time interval between sampling was approximately 2 h. In addition to these dates, samples were collected once a day, with the sampling on November 22, 2021, performed after rainfall. Each sampling involved collecting four L of river water (1 L × 4 replicates, total of four L) for filtration and another 0.25 L for the measurement of basic environmental parameters, i.e., pH, electrical conductivity (EC), and turbidity. After securely tightening the bottle caps, the samples were placed into 45 L transparent plastic bags and then stored in a cooler box to isolate them from samples collected at other locations. Water temperature was measured at sampling sites using a bar-shaped mercury thermometer, and the other parameters were measured from 0.25 L samples transported to the laboratory. EC and pH were measured using a portable multiwater quality meter (HORIBA TOA DKK, Kyoto, Japan) and turbidity was measured using a turbidity meter (PT-200_TROAM, Asahi Science, Tokyo, Japan) in the laboratory. Turbidity, EC, pH, and temperature were not measured in September 2021; temperature was not measured in November 2021; and EC was not measured in June 2023 due to technical problems. During the surveys, tethered buckets were used to collect water samples in cases where entering the river was physically impossible. Each sampling event was completed within 2 h. To prevent sample-to-sample contamination and contamination from external sources during the survey and experiments, a one L of sterile distilled water prepared in the laboratory was placed in a cooler box and transported to the sampling site. After being opened and closed at the sampling site, it was returned to the laboratory with river water samples. We obtained rainfall data for the 24-h period before each sampling event under rainfall conditions from the nearest meteorological station (Akae Observatory, https://www.data.jma.go.jp/obd/stats/etrn/index), and the data are shown in Appendix 1. As a negative control, one L of surface water was collected following a period without rainfall at the same survey sites on December 27, 2021, and November 28, 2022. In addition to one L of surface water, riverbed sediment was collected on August 8 and December 13, 2024.

Filtration and extraction

Figure 2 shows the experimental process. The collected samples were immediately transferred to the laboratory in a cooler box with ice packs, and filtration was completed within 5 h using a glass fiber filter with a pore size of 0.7 µm (GF/F; GE Healthcare Japan, Tokyo; referred to as “small filter”). To avoid filter paper clogging when filtering turbid water, we also used a filter membrane with a pore size of 2.7 µm (GF/D; GE Healthcare Japan, Tokyo; referred to as “large filter”), which is the largest pore size among the glass fiber filters. We compared the filtration efficiency and eDNA recovery rate from each filter with two replicates per filter. The samples collected without rainfall were filtered using only GF/F (one replicate). The funnels and bases (Advantec, Taipei, Taiwan) used for filtration were sterilized in advance via autoclaving (LSX-700, TOMY). Before processing the samples, the funnels, bases, clamps, and tweezers were sterilized again using 10% bleach. We shook the bottle before filtering to prevent eDNA from settling at the bottom due to prolonged filtration time. When the filtration rate dropped below 20 drops/10 s, the filter was replaced with a new one, with a maximum of two filters used per sample. Any remaining unfiltered water was discarded. This limitation was intended to minimize filter clogging and reduce the labor-intensive nature of the extraction process, thereby improving the efficiency of eDNA analysis in turbid water environments (Hinlo et al., 2017). The water volume and filtration time until clogging were recorded for each filter paper. In cases where two filter papers were used for a sample, the total volume and time were used as the volume and time respectively. After processing the samples, all filter papers were frozen at −20 °C until DNA extraction.

DNA extraction was performed using the DNeasy PowerSoil Kit (Qiagen, Hilden, Germany), which produces DNA extracts free of inhibitors (Eichmiller, Miller & Sorensen, 2016). This kit is also able to capture extracellular DNA resulting from unfreezing processes (Gielings et al., 2021; Hermans, Buckley & Lear, 2018). Firstly, using sterilized scissors and tweezers, the filter papers were cut to pieces of approximately 1  × 3 mm on a clean bench that had been pre-treated with ultraviolet sterilization. The workbench was wiped with 70% ethanol between samples to prevent cross-contamination. Subsequently, DNA was extracted from the filter fragments following the kit manufacturer’s protocol and eluted finally in 100-µL volumes. After extraction, the DNA template solutions were stored at −80 °C for subsequent PCR amplification.

We adopted two methods for eDNA extraction from sediment samples. In the first approach, approximately 28 g of sediment sample was mixed with 100 mL of sterile phosphate buffer solution (PBS, pH = 7) for 2 min in a plastic bottle pre-sterilized with bleach for 30 min. Subsequently, 100 mL of the mixture was poured into a sterile funnel, taking care not to introduce sand particles, and filtered through a GF/F glass fiber filter paper (Nevers et al., 2020). In the second approach, 20 g of sediment sample was mixed with 20 mL of phosphate buffer for 15 min. Two milliliters of the mixture was centrifuged (10 min at 10,000 g), and 600 µL of the resulting supernatant was used for subsequent extraction (Rota et al., 2020; Taberlet et al., 2012b). The samples collected in Aug. 2024, and Dec. 2024 were extracted following the DNeasy PowerSoil Pro Kit manufacturer’s protocol and eluted in 100-µL volumes (since the DNeasy PowerSoil Kit has been discontinued).

eDNA quantification and PCR conditions

The species-specific cattle DNA concentration was quantified using primers and probes targeting the mitochondrial cytochrome b (Cytb) region (Dooley et al., 2004). Amplicon size is 116 bp. The primer and probe sequences were as follows: forward primer: CGGAGTAATCCTTCTGCTCACAGT, reverse primer: GGATTGCTGATAAGAGGTTGGTG, Probe: TGAGGACAAATATCATCATTCTGAGGAGCWARGTYA. The mixture was dispensed into the independent wells of a QuantStudio 3D dPCR 20-K chip using a QuantStudio 3D dPCR Chip Loader (Applied Biosystems, Foster City, CA, USA). The end-point PCR was performed using a thermal cycler (ProFlex, Applied Biosystems, Foster City, CA, USA). To determine the optimal PCR conditions, preliminary experiments were performed using a dilution series of DNA extracted from beef meat (10⁻²–10⁻³ fold, starting concentration: 4,738 copies/µL). DNA was extracted from the beef samples using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. The results obtained for the optimal PCR conditions are shown in Appendix 2. When using forward and reverse primers at 450 nM and the probe at 125 nM, an annealing temperature of 60 °C and 50 cycles yielded an adequate separation of fluorescence between positive and negative wells in the dPCR system for both 100-fold and 1,000-fold dilutions. Therefore, all samples were adopted under these conditions for the subsequent PCR amplification. The default reaction volume for QuantStudio 3D dPCR is 14.5 µL, consisting of 7.25 µL 1 × Quant Studio 3D Digital PCR Master Mix, two µL of DNA template solution, 450 nM each primer (the forward and reverse primer) and 125 nM Taqman probe. The fluorescence intensity thresholds defined in the default settings of the Quant Studio 3D Analysis Suite (Applied Biosystems, Foster City, CA, USA) were used to distinguish between positive and negative wells. However, in cases where the fluorescence intensity of positive wells was low and not clearly separated from that of negative wells in the plots, those wells were considered negative (Nukazawa, Akahoshi & Suzuki, 2020). Because the QuantStudio 3D dPCR 20-K chip has been discontinued, non-rainfall samples collected in August and December 2024 were quantified using the QIAcuity one digital system (Qiagen, Hilden, Germany), with a default reaction volume of 40 µL. The same PCR conditions as those mentioned above were used. Since the initial attempts with QIAcuity One dPCR resulted in poor filling of the Nanoplate 26K 24-well plate, the reaction volume was increased to 43 µL following the manufacturer’s recommendation, consisting of 11.25 µL Probe PCR Master Mix, two µL of DNA template solution, 450 nM each primer (the forward and reverse primer) and 125 nM Taqman probe. The mixture was dispensed into a Nanoplate 26K 24-well. Finally, the results were confirmed using the QIAcuity Software Suite. Because the Quant Studio 3D Digital PCR Master Mix and Probe PCR Master Mix may contain bovine serum components, we included two samples without DNA in each type of dPCR to check for reagent contamination (PCR control).

Environmental factors

In this study, specific environmental variables influencing the detectability of eDNA were selected, and their relationship with eDNA concentration was assessed. In the terrestrial environment, the percentages of forest, farmland, and urban areas were taken into consideration, given their relevance to the occurrence of surface runoff (Guzha et al., 2018). Furthermore, two specific distances were considered: the distance from eDNA release sources (the cattle sheds) to the study river (the shortest straight-line, referred to as land distance) and the path that terrestrial eDNA travels to reach the sampling site after being transported into the water body (nonstraight line, referred to as the waterways distance). The land and waterways distances were approximately 125 and 360 m for the Mizunashi River, 36 and 180 m for the Oka River, approximately 311 and 187 m for the Tagami River, and 17 and 134 m for its irrigation channel, respectively (Google Maps). In the aquatic environment, we considered water temperature, pH, EC, and turbidity.

The magnitude and duration of rainfall are potentially key factors regulating terrestrial eDNA runoff. To determine their influence, we fitted the rainfall time series (i.e., hyetograph) for each sampling campaign using a normally distributed probability distribution to describe the characterization of each rain event. The hourly precipitation recorded over a 24-h period at the Akae rainfall station was used as an observation value, and a function describing a calculated precipitation was developed; αf(x), where f(x) and α are the normal distribution function and a constant, respectively. Here, the time when the water sampling started was set to time 24, and 24 h before sampling started was set to time 1 to initialize the time axis; thereby, the mean (µ) of the normal distribution varied from 1 to 24. The search range for the standard deviation (σ) of the normal distribution and the constant α were set from 0 to 100 and from 0 to 1,000, respectively, because these parameters varied greatly in the preliminary attempts. Iterative parameter searches were performed using R ver. 4.3.2 (R Core Team, 2023) to determine the values of μ, σ, and α that minimized the root mean square error of observation and calculated value. A larger μ suggested that the peak rainfall occurred closer to the sampling time, whereas a larger σ indicated that similar magnitude rainfall events lasted longer (Fig. 3, Appendix 1).

Regression analysis

Generalized linear mixed models (GLMM) with the Gaussian family were used to evaluate the relationship between eDNA concentrations and environmental factors (μ, σ, turbidity, land distance, waterways distance, percentage of each land use, pH, EC, water temperature, and type of filter), with each of the four rivers as a random effect. In the models, the interaction between turbidity and the type of filter paper was considered based on a previous report that examined the effect of suspended solids on eDNA particle distribution size (Barnes et al., 2021). Prior to modeling, the collinearity between environmental factors was determined using generalized variance inflation factors (GVIF) using the R package “car” (Fox & Weisberg, 2019). The factor with the highest adjustment value of GVIFˆ(1/(2*Df)) was repeatedly omitted until all fell below 2 (Fox & Monette, 1992). Through this procedure, the percentages of the forest, farmland, and urban areas and water temperature were excluded from the model. The “glmmTMB” package (Brooks et al., 2017) was used to build the GLMMs using log-transformed eDNA concentration as the dependent variable to identify factors affecting eDNA content. One was added to eDNA concentration prior to log-transformation to accommodate values of 0 (Everts et al., 2024). All statistical analyses were conducted in R ver. 4.3.2 (R Core Team, 2023). Additionally, two generalized linear models (GLMs, gaussian family) were used to evaluate the relationships between filtration time and turbidity and filter paper type, and the relationships between filtration volume, turbidity, and filter paper type. To determine whether the concentration of eDNA differed between small filter and larger filter, we used the paired Wilcoxon rank-sum test.

Measurements of basic water quality parameters

Table 1 shows the measurements of the basic water quality parameters in each river. The ranges of water temperature at each river are as follows: Mizunashi, 14.90 °C–18.00 °C (average, 17.30 ± 0.60); Oka, 18.40 °C–21.60 °C (average, 20.32 ± 0.57); Tagami, 19.00 °C–24.20 °C (average, 21.30 ± 1.08). The water temperature of irrigation canal was recorded only once in November 2022, and it was 19.00 °C. The pH values at each river are: Mizunashi, 6.61–7.96 (average, 7.54 ± 0.48); Oka, 6.09–7.70 (average, 7.18 ± 0.63); Tagami, 5.98–7.27 (average, 6.95 ± 0.49); irrigation canal, 7.08–7.33 (average, 7.20 ± 0.12). Based on observations performed twice daily, the pH values of the samples at the same site tended to decrease over time. The EC values at each river are as follows: Mizunashi, 53.60–119.30 μS/cm (average, 83.40 ± 29.69); Oka, 55.70–76.70 μS/cm (average, 67.70 ± 8.91); Tagami, 39.80–126.60 μS/cm (average, 85.96 ± 36.52); irrigation canal, 25.80–166.20 μS/cm (average, 118.90 ± 80.63). The fluctuations in EC in Tagami and irrigation canal were larger than those in Oka and Mizunashi. The turbidity at each river is as follows: Mizunashi, 5.78–285.82 ppm (average, 78.13 ± 97.69); Oka, 46.73–314.34 ppm (average, 142.23 ± 108.30); Tagami, 4.07–138.54 ppm (average, 31.04 ± 52.86); irrigation canal, 0.45–189.94 ppm (average, 63.78 ± 109.26). Turbidity fluctuated greatly in all sites.

Filtration volume and time

The filtration volume and time for each DNA concentration method are shown in Appendix 3. A higher filtration volume and a shorter filtration time were observed for the large filter than for the small filter. For the samples collected from November 2021 to June 2023, two filter papers were used in 87.71% of the small filter-based filtrations. Among them, the minimum turbidity detected was 4.07 ppm. Additionally, for samples with turbidity equal to or greater than 25.38 ppm, even two filter papers did not allow full filtration (61.10% of the samples). The sample collected during the survey conducted in November 2022 had the longest filtration time (filtration volume: 690 mL, filtration time: 2,610 s). For large filter-based filtration, two filter papers were used for 42.86% of the samples. For samples with turbidity equal to or greater than 110 ppm, a water volume of one L could not be filtered even using two filter papers (72.22% of the samples). Overall, large filters were more efficient in filtering high-turbidity samples than small. Table 2 shows the relationships between filtration volume, filtration time, filter type, and turbidity. The GLMs indicated that turbidity had a significant positive effect on filtration time, and that filter type had a significantly greater effect on increasing filtration time. In addition, turbidity had a significant negative effect on filtration volume, and filter had a significantly negative effect on filtration volume.

eDNA quantification

Target DNA was not detected in the autoclaved distilled water sample transported to the site on each survey date and the PCR control, confirming the absence of contamination between samples during the surveys and reagent contamination. Figure 4 shows the results of the eDNA survey conducted in the rainfall and non-rainfall events. In non-rainfall events, the target eDNA was not detected in most samples, and in the few where it was detected, the concentration was extremely low (0–338 copies/L). In contrast, considerably more target DNA was observed during the rainfall period (November 2022 samples).

Figure 4 shows the results of eDNA quantification during the rainfall period. Environmental DNA was detected in 42 out of 47 samples. In terms of survey periods, the eDNA concentrations were in the range of 855–8,267 copies/L for September 2021, 0–1,350 copies/L for November 2021, 250–3,801 copies/L for May 2022, 0–2,730 copies/L for September 2022, 4,807–399,185 copies/L for November 2022, and 0–1,040 copies/L for June 2023. In addition, during the heavy rainfall event of November 29, 2022, high concentrations of cattle DNA were detected, whereas in November 2021, which was a less rainy month, many samples showed low eDNA concentrations or even none. Among the five samples with no eDNA traces collected throughout the entire investigation period, three were from the Mizunashi River. During the survey period on the same date, higher eDNA concentrations tended to be detected in the irrigation channel than in the Tagami River.

Relationships between environmental factors and eDNA concentration

The GLMMs suggested significant positive effects of σ, turbidity, pH, and filter and significant negative effects of waterways distance on eDNA concentration. However, no significant effects of μ, land distance, EC, or the interaction between turbidity and filter were observed (Table 3).