Impacts of Flood Disturbance on the Dynamics of Basin-Scale Swimming Fish Migration in Mountainous Streams

Abstract

This study investigated the response of sweetfish (Plecoglossus altivelis), a species that migrates dynamically throughout a river basin, to two flood events during the summer of 2020 in the Nagara River, located in the central region of Japan. By combining multiple environmental DNA (eDNA) surveys and hydrological modeling, the spatiotemporal distribution of P. altivelis throughout a mountainous river basin was captured and analyzed. The eDNA concentrations at 42 sites in the Nagara River Basin were analyzed five times from August to early October 2020. In addition, Rainfall–Runoff–Inundation model calculations were performed using 1-km resolution precipitation data as input values to analyze the magnitude of the flood disturbance at the eDNA sites. The daily specific discharge Qs (m³/s/km²) was employed as an index of the flood magnitude. The calculation period included Flood Events 1 and 2 with 52 and 38 days of precipitation reaching 1923 and 528 mm, respectively. The results of the eDNA analysis showed that, immediately after Flood Event 1, the spatial distribution of P. altivelis was unevenly distributed in the upper reaches of the Nagara River and some of its tributaries. Subsequently, the distribution expanded to the entire mainstem. The distribution of the maximum daily specific discharge suggested that the river segments with high eDNA concentrations of P. altivelis immediately after the high-magnitude flood event were those with a relatively low intensity of flood disturbance compared with those in the other connected river segments and tributaries. The results of this study indicate that the resilience of riverine communities to extreme floods is supported by the continuity and connectivity between the mainstem and its tributaries in mountainous river basins.

1. Introduction

Flow fluctuation is one of the most universal phenomena in rivers [1] and is a critical factor that determines the characteristics of river ecosystems [2]. Human activities, such as dam construction, affect natural flow regimes, with their effects varying depending on the species of plants and animals [3]. Flood disturbances have a strong impact on the diversity of riverine ecosystems. From mountain streams to lowland rivers, riverine fishes have a life history adapted to flooding [4,5]; however, extreme flood events can also eliminate individual organisms and significantly alter the biota [3,6,7,8]. Riverine fishes are known to use various locations as refuges during floods [9]. In a river with a gentle gradient, the inundation area with a low flow velocity is used as a refuge [10], whereas, in a steep mountainous stream, the external forces acting on the streambed by floods are large, limiting the space for organisms to take refuge. Stable boulders in streams may serve as refuges from floods [11]. In addition, tributaries serve as important refuge areas during mainstream floods [12]. Although some information at limited spatial scales is available on flood refuges for fish in mountainous streams, the migration patterns of highly mobile fishes within a river basin in response to flood disturbance have not been investigated at the river basin scale.

In recent years, the environmental DNA (eDNA) analysis method has made it easier to determine the distribution of species using simple water sampling [13,14]. In addition, a quantitative evaluation of the abundance of freshwater fishes using eDNA analysis has been reported [15]. However, eDNA concentration is an unstable indicator that fluctuates depending on the environmental conditions [16]. To properly interpret the results of the eDNA analysis, the complex and dynamic environments from which eDNA is sampled need to be considered for accurate biological interpretation [17]. One solution to this problem is the integration of eDNA analysis with hydrological observations and models, which has the potential to create new values for understanding and better managing river ecosystems [17,18].

Therefore, we proposed a pioneering approach that combines multiple eDNA surveys and hydrological modeling to capture and analyze fish migration in response to flood disturbance throughout a mountainous river basin, focusing on the sweetfish Plecoglossus altivelis as a model species. Multiple eDNA surveys were conducted before and after flood events to capture the response of organisms to flood disturbances. To accurately estimate the degree of flood disturbance acting on each river segment in a river system, hydrological model calculations were performed using high spatiotemporal resolution precipitation data as the input values. This combination was applied to determine the relationship between flood disturbance and the spatiotemporal distribution of highly mobile organisms at the river basin scale, illustrating how the swimming fish P. altivelis responded to the flood event that occurred during the summer of 2020 in the Nagara River Basin, located in the central region of Japan. The results of this study will help to enhance our understanding of the resilience of riverine communities to floods. In addition, this approach may be applicable for analyzing the impact of climate change and human activities on the resilience of river ecosystems.

2. Materials and Methods

2.1. Study River

The study was conducted in a middle–upper segment > 42.0 km from the mouth of the Nagara River (mainstem) and in its four main tributaries (Yoshida, Itadori, Mugi, and Tsubo Rivers) in Central Japan (Figure 1). The entire channel length of the Nagara River is 166 km, and its basin area is 1985 km². The highest point within the watershed is Mt. Dainichi at 1709 m. The longitudinal profiles of the mainstem and the four major tributaries are shown in Figure 2. The sections at approximately 42.0–53.0, 53.0–75.0, and >75.0 km of the Nagara River flow through the alluvial fan, valley bottom plain, and mountain area, respectively. The average bed slope of the alluvial fan and valley bottom plain section is approximately 1/560. The bed slope of the mountain area section ranges from 1/250 to 1/10. Among the four main tributaries, the Yoshida River and the upstream sections of the Itadori River have steeper slopes. The bed material in most of the sections where eDNA surveys were conducted consisted of boulders and gravels. In the section from the valley bottom plain to the alluvial fan of the mainstem, gravel–sand bars are formed and widely appear under normal water conditions. Boulders are particularly common in the upper reaches. There were no dams in the study section of the mainstem of the Nagara River. Therefore, upstream fish migration was ensured through the mainstem. Several check dams with fish passages were present in the Yoshida and Mugi Rivers.

This region and P. altivelis were designated a “globally important agricultural heritage system” by the UN Food and Agricultural Organization in 2015 for their value in combining livelihood, nature, and culture [19]. Plecoglossus altivelis is an amphidromous fish that uses rivers and coastal areas during its annual lifespan. They spend their larval to juvenile stages in estuaries and coastal areas in the winter and migrate to upstream rivers in the spring [20,21]. Furthermore, they have one of the highest swimming abilities among swimming fish. The ratio of their critical swimming speed is 3.1–8.0 times their body length, and the mean ratio is approximately five times their body length [22]. These fish use the middle and upper reaches of rivers as their main habitat during summer growth and feed on algae attached to boulders and cobbles. Mid-October to the end of November is the spawning season in the recent Nagara River, and the main spawning site is located in the alluvial fan section. There are six fishery cooperatives in the study rivers. They release plenty of hatchery-raised P. altivelis juveniles (approximately 10 cm in standard length) in each of their territories every year from April to June, mainly in May. Therefore, both natural and hatchery fish are widely distributed in the Nagara River system during the early summer.

2.2. Flood Events

Japan is located in the Asian monsoon region, where relatively high flows are observed between June and July (East Asian rainy season) and between September and October (typhoon season). From November to May, the flows are relatively low and stable. The precipitation in 2020 was record-breaking, especially in June and July, when torrential rainfall was observed over a wide area of Japan.

The precipitation observation record at the Hachiman rain gauge located upstream of the Nagara River Basin (Figure 1) is shown in Figure 3. In the 52 days from 10 June to 30 July, 1923 mm of precipitation was observed, which was more than three times that of the normal June and July precipitation at the site. In the Nagara River, intermittent rains during this period caused flooding that lasted for approximately two months. The Yoshida River, one of the main tributaries, experienced heavy runoff in early July (Figure 4). Precipitation was low in August, and strong rainfall due to typhoons was observed several times in September. In total, 528 mm of precipitation was observed in 38 days, from 2 September to 10 October. In this study, these two flood events were defined as Flood Events 1 and 2, respectively.

The effects of “extreme” Flood Event 1 and less severe Flood Event 2 on P. altivelis, which were distributed throughout the Nagara River Basin before the flood event, were analyzed by eDNA survey and a hydrological model analysis.

2.3. eDNA Sampling

The sampling sites were chosen to examine the spatiotemporal distributions of P. altivelis in the Nagara River system immediately after Flood Event 1. Water was sampled for the eDNA assay from 42 sites (mainstem, 20; Yoshida River, 6; Itadori River, 3; Tsubo River, 5; and Mugi River, 8) on each of five sampling days: 11 August, 25 August, 14 September, 23 September, and 6 October (Figure 1). It took two days (those above and the following days) to sample water from the 42 sites. The sampling period corresponded to the summer growing season and just before the autumn spawning season.

We collected 1-L water samples using a polypropylene bottle from the water surfaces at a location with a current velocity of approximately 50 cm/s. Benzalkonium chloride solution (0.1% w/v) was added to the bottles to preserve the eDNA [23]. The samples were stored in a cooler box, along with a 1-L distilled water blank, placed in a laboratory, and refrigerated.

2.4. DNA Filtration, Extraction, and Quantitative Real-Time PCR

The water samples were pulled by vacuum through a 0.7-µm membrane filter (47-mm GF/F glass fiber; GE Healthcare, Little Chalfont, UK) within 2 days of sampling. The cooler blank was then filtered simultaneously. A filtered “equipment blank” (1-L DNA-free water) was incorporated before and after each sample filtration day. The “cooler blank” and “equipment blank” were negative controls, and their DNA levels were measured to identify any DNA contamination from field preparation and transportation, filter equipment, or background sources. The filters were placed in a light-shielding bag and stored at −20 °C until eDNA extraction. DNA extraction was performed according to the method described by Doi et al. [14]. The final volume of each extracted DNA sample was 100 μL, and all samples were stored at −20 °C until the PCR assay.

The eDNA of P. altivelis was quantified via real-time TaqMan PCR (qPCR) using the PikoReal Real-Time PCR system (Thermo Fisher Scientific, Waltham, MA, USA). Mitochondrial cytochrome b gene fragments (131 bp) were amplified and quantified with a TaqMan probe using the following primers: Paa-CytB-Forward: 5′-CCTAGTCTCCCTGGCTTTATTCTCT-3′, Paa-CytB-Reverse: 5′-GTAGAATGGCGTAGGCGAAAA-3′, and Paa-CytB-Probe: 5′-FAM-ACTTCACGGCAGCCAACCCCC-TAMRA-3′. The species-specificity of the primer/probe set was confirmed by Yamanaka and Minamoto [13]. Each TaqMan reaction mixture contained 900 nM of each primer (Paa-CytB—Forward and Reverse), 125-nM Paa-CytB-Probe, 10-µL qPCR Master Mix (TaqMan Environmental Master Mix 2.0; Life Technologies, Carlsbad, CA, USA), 0.4-µL uracil N-glycosylase (Thermo Fisher Scientific), and 4-µL DNA solution. The total reaction mixture volume was 20 µL, which included DNA-free water. The PCR conditions were as follows: 50 °C for 2 min, 95 °C for 10 min, and 55 cycles of 95 °C for 15 s and 60 °C for 1 min. There were three replicates per water sample. Four no-template controls were included in each PCR run. The limit of detection of the qPCR was one copy per reaction, with an average of three replicates. A standard curve was constructed using a dilution series of 10,000 (n = 2), 1000 (n = 2), 100 (n = 2), and 10 copies (n = 2) per PCR reaction. The R² values of the standard curves ranged from 0.956 to 0.995. The standard samples were the target DNA cloned into the plasmids.

2.5. Rainfall–Runoff–Inundation Model

The Rainfall–Runoff–Inundation (RRI) model [24] was employed to calculate the flow discharge time series at any point on a river that constitutes the Nagara River system. The model is a two-dimensional grid-based model that considers water movement on both slopes and river channels. The model calculation requires rasterized topographic information, as well as distributional rainfall input. In this study, HydroSHEDS [25] with 15-s (approximately 450 m) resolution, provided by the World Wildlife Fund [26], was applied for topographic information. Radar/Rain gauge-Analyzed Precipitation (RA) data from the Japan Meteorological Agency were used as the rainfall input. This dataset is created by combining radar precipitation and ground precipitation, and it represents 1-km grid precipitation with some accuracy [27], even in areas where rain gauges have not been installed. The RRI model, a type of distributed runoff analysis model, was used for the calculations in this study. Any model that can calculate the river discharge at any point in the watershed using the spatiotemporal distribution of rainfall as an input value can be adopted as an alternative method.

The RRI model requires the operator to set the parameters for the calculations. Most of the variables are provided by RRI_BUILDER, a GUI application for setting up RRI models. In this study, nine sets of parameters were set for each grid, depending on the combination of land cover and soil, which were automatically set by RRI_BUILDER. Calibration was performed using RA data from 2018 to set the appropriate parameters. The model parameters were based on the values given in the RRI model manual [24], and only the Manning’s roughness on the slope cells and soil depths were varied in several patterns to find the preferred parameter settings (

Table 1. List of parameters of the Rainfall–Runoff–Inundation model.

Explanation	Value	Determination
Manning’s roughness in river channel	0.030	Default value
Manning’s roughness on slope cells	0.500	By calibration
Soil depths (m)	1.000	By calibration

).

The calculation results were output as a flow time series at a specified point. The calculation results were calibrated and validated by calculating the Nash–Sutcliffe model efficiency coefficient (NSE) and the coefficient of determination (R²) using the daily discharge data from the Japan Ministry of Land, Infrastructure, Transport and Tourism at the Chusetsu water gauging station, located near the downstream of the study area (Figure 1). The calculation period was from January to December 2020. However, because the RRI model cannot handle both snow accumulation and snowmelt, only the calculation results from April to November, excluding the period when snow accumulation was observed in the Nagara River Basin, were considered in this study.

2.6. Hydrologic Index of Flood Magnitude

The disturbance regime caused by flow fluctuation can be treated quantitatively by decomposing it into five components: magnitude, frequency, duration, timing, and rate of change [1,28]. Olden and Poff [29] calculated 171 hydrological indicators from 420 sites and 36 years of flow data in the U.S. and identified 26 important indicators. These hydrological indices are calculated from the daily discharge (m³/s). As the flow discharge of a river is affected by the size of the river, there are methods to ensure that the size of the river does not affect the analysis results, such as dividing the flow rate by the average value or dividing the flow rate by the catchment area and treating it as a specific discharge [29]. As the purpose of this study was to examine flood disturbances at several sites within the same river basin for a relatively short period of time, specific discharge was used to calculate the hydrological indices in this study.

Among the five components of flood disturbance, we focused only on the magnitude. The reason for this was to examine the effect of the flood event immediately before the survey on the spatial distribution of the eDNA concentration. The maximum daily specific discharge Q_smax (m³/s/km²) during the flood event was used as an indicator of the flood magnitude.

2.7. Analysis

The relationships between the eDNA concentrations at 42 eDNA sites and the index of flood magnitude were examined for two flood events. The eDNA surveys conducted on 11 August and 6 October were used to determine the spatial distribution of P. altivelis after the two flood events. As the variation in eDNA concentration had a range of several orders of magnitude, it was treated as an index on a logarithmic scale. To understand the impact of flood disturbance on the spatial distribution of P. altivelis, we aggregated the 42 eDNA sites into eight segments; that is, four mainstream segments of the Nagara River (S1–4, numbered from upstream) and four main tributaries (Yoshida, Itadori, Mugi, and Tsubo Rivers). The relationship between eDNA concentrations and flood disturbance indices in the eight river segments was examined.

3. Results

3.1. Spatiotemporal Change in eDNA Concentrations of P. altivelis

The eDNA concentrations of P. altivelis changed spatially and temporally in the Nagara River Basin (Figure 5). On 11 August, after Flood Event 1, relatively high eDNA concentrations were detected in the upper segment of the mainstem and in the Itadori River. On 25 August, relatively high eDNA concentrations extended to the mainstem immediately downstream of the confluences with the Yoshida and Itadori Rivers.

During September, which was mostly included in Flood Event 2, eDNA tended to be distributed throughout the mainstem, although the eDNA concentrations were particularly low due to dilution caused by high river flows. On 6 October, the eDNA concentrations were evenly high throughout the mainstem in normal flow after Flood Event 2.

3.2. RRI Model Calculation Results

The river flow of the Nagara River system in 2020 was calculated with high accuracy using the RRI model with RA data. The daily discharge hydrograph observed at the Chusetsu water gauging station and the RRI model calculation results at the same location are shown in Figure 6. The calculated discharge was slightly underestimated for Flood Event 1 with very heavy rainfall. Conversely, for Flood Event 2, which was a small-to-medium-sized flood event, the calculated results were slightly larger than the observed values. The NSE and R² values were 0.83 and 0.88, respectively, indicating that the reproduced calculation results were sufficiently accurate. The spatial resolution of the RA data used as input for the RRI model calculation was 1 km, suggesting that the flow time series of each tributary constituting the Nagara River system was also calculated with sufficient accuracy.

3.3. Flood Magnitude Index

The time series of the daily discharge at the 42 eDNA sites calculated by the RRI model was divided by the watershed area A (km²) of each site and converted into a time series of daily specific discharge Q_s (m³/s/km²). The analysis periods of the flood magnitude index corresponding to the two flood events were 62 days from 10 June to 10 August for Flood Event 1 and 56 days from 11 August to 5 October for Flood Event 2 (Figure 6). The maximum daily specific discharge Q_smax during the flood events was calculated.

The Q_smax values for the two flood events are shown in Figure 7. The maximum value of Q_smax in Flood Event 1 was 2.11 at the Yoshida River, and the mean value for the 42 sites was 1.37 (Figure 7a). Focusing on the four main tributaries, it was confirmed that the Yoshida River had the highest flood magnitude, followed by the Tsubo River, and the Itadori and Mugi Rivers had relatively low magnitudes. On the mainstem segments of the Nagara River, it was confirmed that the flood magnitude was higher downstream of the confluence with the Yoshida River than that upstream of the confluence. The maximum value of Q_smax for Flood Event 2 was 0.75 at the Mugi River, and the mean value of the 42 sites was 0.50 (Figure 7b). The flood magnitude of the mainstem segment of the Nagara River tended to be slightly lower than that of the main tributaries. The overall flood magnitude was lower than that of Flood Event 1.

3.4. Relationship between eDNA Concentration and Flood Magnitude Index

The relationships between the eDNA concentrations and average Q_smax in the eight river segments examined are shown in Figure 8. The mainstem segments of the Nagara River were numbered 1–4 from upstream to downstream. The connection between the four main tributaries and the four mainstem segments is also shown in the schematic image in Figure 9.

The mean eDNA concentration varied among the river segments after Flood Event 1 (Figure 8a). The relationship between Nagara S1 and the Yoshida River located in the upper reaches and Nagara S2 located downstream of their confluence illustrated that the higher the flood magnitude, the lower the average eDNA concentration. Similarly, focusing on the relationship between Nagara S2 and the Itadori River and Nagara S3 downstream of their confluence, higher average eDNA concentrations were found in the Itadori River, which had the lowest flood magnitude. After Flood Event 2 (Figure 8b), a different trend was observed. The eDNA concentrations in the entire mainstem of the Nagara River and in the Itadori River became uniform. Conversely, the eDNA concentrations in the Yoshida, Mugi, and Tsubo Rivers, which were low immediately after Flood Event 1 (Figure 8a), remained low after Flood Event 2.

4. Discussion

The combined approach of multiple eDNA surveys and hydrological modeling showed that it is possible to capture and analyze fish migration in response to flood disturbances in the mountainous river basin. The biased spatial distribution of P. altivelis after Flood Event 1 with record precipitation could be explained by the magnitude of the flood disturbance. This result suggests that the fish moved between the directly connecting segments to evacuate from prolonged and high-intensity floods (Figure 8a). After that, P. altivelis were widely redistributed in the segments of the mainstem, presumably because the mainstem has more available habitat than the tributaries. Moreover, floods of lesser magnitude might not affect the spatial distribution of P. altivelis (Figure 8b). In addition, several small tributaries flow into the mainstem of the Nagara River. It is possible that individuals who took refuge in these tributaries returned to the mainstem, where they could find more habitats. The behavior of the P. altivelis population evacuated from the swollen segment to the other segment was compatible with the observations by Koizumi et al. [12]. However, there was insufficient material to determine whether floods induce migratory behavior in P. altivelis, using the magnitude of the maximum specific discharge as a threshold. Although the specific discharge is useful for making comparisons between sites with different catchment areas, the magnitude of the external forces acting on the streambed cannot be defined solely by the specific discharge.

The response to floods may differ among fish species based on their ecology and biological characteristics. Plecoglossus altivelis is known to avoid turbid water and evacuate turbid sections [30]. The Yoshida River, which had the lowest eDNA concentration after Flood Event 1, showed high turbidity during the flood (Figure 4). In addition, turbidity inhibits the growth of attached algae, which is the dominant food resource for P. altivelis. Therefore, if the flood period is long, P. altivelis populations may move as a “habitat shift” due to a decline in attached algae rather than “evacuation” from floods. Thus, some ecological factors may be contained in our results of the relationship between flood disturbance indices and eDNA concentrations. Future studies integrating the approaches of this study and conventional ecological works can highly progress the understanding of riverine fish ecology.

The results of this study indicate that the resilience of riverine communities to extreme floods is supported by the continuity and connectivity between the mainstem and its tributaries in mountainous river basins. The resilience of riverine ecosystems will become more important as the magnitude and frequency of floods are projected to increase as a result of climate change [31]. Hydrological model calculations using climate model products as input values may reveal the impact of climate change on riverine ecosystems via its impact on flow regimes. The integration of eDNA analysis with hydrological observations and models has the potential to create new values for understanding and better managing river ecosystems.