Distribution of Japanese Eel Anguilla japonica Revealed by Environmental DNA

The abundance of Japanese eel Anguilla japonica has rapidly decreased in recent decades. Following a re-evaluation of the possibility of extinction, the Japanese Ministry of the Environment and the International Union for Conservation of Nature listed the Japanese eel as an endangered species in 2013 and 2014, respectively. However, their abundance and precise distribution have never been clarified owing to their nocturnality and difficulty in their capture. In this study, the distribution of Japanese eels was investigated by monitoring for environmental DNA (eDNA), a non-invasive and efficient detection method. A total of 365 water samples were collected from 265 rivers located throughout Japan. High concentrations of eDNA of Japanese eels were detected in rivers on the Pacific side, but were low in the Sea of Japan side. In particular, very little eDNA amplification was confirmed from Hokkaido and the north of the Sea of Japan. The eDNA distribution in Japanese rivers coincides with the transport of the larvae in the ocean, as estimated by numerical simulations. Generalized linear mixed models were developed to explain the distribution of eDNA concentrations. The total nitrogen concentration emerged as an important factor in the best model. These results indicate that the distribution of Japanese eel is mostly determined by the maritime larval transport, and their survival and growth depend on the abundance of food in the river. The findings of the present study are useful for the management of populations and in the conservation of Japanese eels.


INTRODUCTION
The Japanese eel (Anguilla japonica) is a catadromous species distributed in the western Pacific Ocean. It is a commercially important fishery species in East Asia and is very popular as traditional seafood in Japan. However, the catch of naturally occurring eels in Japanese rivers has rapidly decreased in recent decades, indicating a remarkable decline in its biomass. In response to this situation, the Japanese eel was classified as an endangered species in Japan by the Ministry of the Environment (2013). The International Union for the Conservation of Nature has also classified it as Endangered in the Red List of Threatened Species (Jacoby and Gollock, 2014). This decline could be related to factors such as overfishing, loss of habitat (Tatsukawa, 2003), and changes in the ocean environment affecting larval survival and recruitment (e.g., Kimura et al., 2001). However, the specific cause remains unclear.
Japanese eels spawn in the west of the Mariana Ridge (e.g., Tsukamoto, 2006). After hatching, the larvae (leptocephali) are transported westward by the North Equatorial Current for more than thousands of kilometers from the spawning site. The North Equatorial Current bifurcates into the northward Kuroshio Current and the southward Mindanao Current at its westernmost boundary off the east coast of the Philippines. Larvae enter the Kuroshio Current in the bifurcation zone, where they metamorphose into glass eels, and are transported to the coast (Figure 1). They then swim toward nearby estuaries and rivers for further growth (Zenimoto et al., 2009;Han et al., 2012;Hsiung et al., 2018). Although it has been suggested that some eels spend their entire lives in seawater without experiencing freshwater environments (Tsukamoto and Arai, 2001), freshwater and brackish waters are generally the main habitats for this species. A recent study on otolith strontium isotope ratios of adults collected from the spawning area has shown that they inhabited Japanese rivers and/or estuaries in the juvenile stage (Otake et al., 2019). This result strongly supports the idea that Japanese rivers and/or estuaries act as nurseries for Japanese eels and are hence important for their reproduction. It is essential to note that the lower reaches of rivers are the main nurseries for eels because they spend the longest part of their life in these habitats, and humans can manage and preserve these environments.
However, much remains unknown about the ecology of the Japanese eel. Despite its acknowledged existence, the life history of eels in the sea is still largely unknown. Research into the mechanism of spawning migration is not yet complete, although some progress has been made. In addition, Japanese eels migrate to East Asian countries, but little research has been done on their detailed distribution within those countries. To promote effective conservation measures within the constraints of human resource, time, and money, it is essential to determine their distribution within specific political boundaries and to identify areas where conservation measures should be prioritized.
One of the major difficulties in determining or confirming the distribution of an aquatic species such as the eel is the low detection rate when using conventional methods. Backpack electrofishers are often used for monitoring eels in freshwater areas (Reid, 2011), but they require advanced skills to collect eels from turbid waters. Eels tend to hide in refuges such as holes and crevices or burrows in the mud during the day (Aoyama et al., 2005), which further contributes to the difficulty in their capture. Furthermore, electrofishers are ineffective in deep areas and are disabled in saline water. Consequently, accurate data on the spatial and temporal distribution of Japanese eels are lacking. Alternative methods are required to monitor their precise distribution and abundance.
Animals and plants release their cellular material including DNA into aquatic or terrestrial environments through excrement, body mucus, blood, and sloughed tissues. This DNA suspended in environmental samples, such as water and sediment is called environmental DNA (eDNA, Miya et al., 2020). By filtering a certain amount of natural water, eDNA is concentrated and captured on a filter, from which it is extracted and subjected to various molecular biology experiments for detection of organisms (Deiner et al., 2015;Turner et al., 2015;Miya et al., 2016). Recent studies have demonstrated that eDNA detection can be a reliable method for determining the distribution of various species of fish in freshwater ecosystems (e.g., Minamoto et al., 2012;Thomsen et al., 2012;Yamanaka and Minamoto, 2016). eDNA detection methods are more sensitive than traditional sampling methods, such as electrofishing or visual surveys, particularly when determining the presence of species that are found at low densities (Jerde et al., 2011). Takeuchi et al. (2019) successfully detected eDNA in the southern west Mariana Ridge, which is one of the most plausible spawning areas of Japanese eels, but spawning eels have rarely been captured there (Kurogi et al., 2011). Itakura et al. (2019) conducted surveys of Japanese eels in several small rivers using both electrofishing and eDNA analysis. In their study, the eDNA of Japanese eels was detected at 56 (91.8%) of the 61 study sites from which individuals were collected by electrofishing, and at an additional 35 sites where individuals were not directly collected. This indicates that eDNA analysis was more sensitive for detecting the presence of eels than electrofishing. They also indicated the possibility of estimating eel abundance and biomass from eDNA concentrations, by showing a weak but significant positive relationship between the eDNA concentration and the abundance and biomass of eels.
It is also difficult to survey many places in a short period using conventional methods. Since eels are an important fishery resource, habitat surveys have been conducted in various rivers by many Prefectures (e.g., Research Institute of Environment, Agriculture and Fisheries, Osaka Prefecture, 2014). However, each survey is limited to less than several rivers in each target region, and there are no studies that cover the entire country. On the other hand, the Ministry of Land, Infrastructure, Transport, and Tourism conducts the National Census on River Environments to monitor the environment and living organisms of 109 first-class rivers across the county. The census provides valuable information on the distribution of organisms, including eels. However, it takes 5 years to survey all the rivers. Nakamura et al. (2015) mentioned that the number of survey points in each river is small, and it is necessary to increase the sampling effort for obtaining detailed data on the environment and biodiversity of the rivers, and based on this data, appropriate river improvements should be undertaken. However, it is difficult to increase the scope and frequency of surveys because of the costs (for labor and to undertake the surveys) involved (Fujita, 2018). Recently, the applicability of the eDNA metabarcoding assay for fish surveys was investigated using case studies of the National Census on River Environments (Kitagawa et al., 2020). The checklist of fish compiled using eDNA metabarcoding and traditional sampling methods were compared, and the results showed good agreement, indicating the effectiveness of the eDNA method for fish monitoring. These recent studies suggest that eDNA could be a powerful tool for monitoring the spatial distribution of eels across large landscapes.
Therefore, in this study, we used the eDNA method to detect the distribution of A. japonica in Japan. We collected water samples from rivers across the country and analyzed them for eDNA. As highlighted above, Japanese eels lay their eggs in the west of the Mariana Ridge, and larvae are transported to East Asia. The transport process and the location of the juveniles transported by the currents would be important in determining the distribution of eels in Japan. Therefore, numerical simulations were used to investigate the transport process of eel larvae to the Japanese coast. Survival after reaching the coast could be dependent on the environment of each river. Therefore, the data on various environmental factors in rivers were compiled using a Geographic Information System (GIS). The results of the eDNA analysis were compared with the results from the simulation, and modeled with the environmental factors to identify the factors determining the distribution of the Japanese eel.

Study Areas
To analyze the eDNA of Japanese eels, samples were collected from the lower reaches of rivers located throughout Japan (Figure 1). Based on climatological and oceanographic insights, the studied area was divided into eight regions: the Nansei-Shoto Islands (NSI), west of Kyushu (WK), Pacific side of the main island of Japan (PMJ), the Seto Inland Sea (SIS), the western part of the Sea of Japan (WSJ), northern part of the Sea of Japan (NSJ), Pacific side of Tohoku (PT), and Hokkaido (HD) ( Figure 1A). NSI is located in the most southern area, and the Kuroshio flows northeastward along these islands ( Figure 1B). This area is the first location where eel larvae reach. However, the rivers in NSI are very short and small because of the narrow area and flat terrain. After passing through the NSI, the Kuroshio mainstream flows eastward in the south of Japan along the Kyushu, Shikoku, and western part of the Honshu area, and moves away from Japan at Cape Inubozaki. In WK, the Tsushima Current diverges from the Kuroshio and flows northward. The Tsushima Current enters the Sea of Japan via the shallow Tsushima Strait (∼120 m depth) and flows eastward along the western Honshu. In NSJ, the Tsushima Current shows unstable flows away from the coast. A weak branch of the Tsushima Current passes from the Sea of Japan to the Pacific Ocean through the Tsugaru Strait and turns south along the coast of the Pacific side of northern Honshu. However, PT is mainly influenced by the cold Oyashio Current. A warm water mass (warm core ring) sometimes separates from the Kuroshio and approaches PT. HD is in the northernmost part of the study area so that the Kuroshio and Tsushima Currents have little effect. It is located in the subarctic zone and water temperatures in the HD rivers are low and severe for eels, especially in winter. SIS is the largest estuary in Japan and connects the Pacific Ocean via the narrow Bungo Channel and Kii Channel. The seawater in SIS is originally from the Pacific Ocean, so that water temperature is moderate, but strongly affected by the terrestrial water.
It was expected that the eDNA of eels would be detected in high concentrations in the PMJ, which is directly under the effect of the Kuroshio Current, and where eel larvae are easily transported. After reaching the estuaries, water quality influences the survival and growth of eels, which are represented by the eDNA concentration in each river.

eDNA Field Sampling
A total of 365 study sites were selected from the downstream reaches of 265 rivers spatially spread across Japan ( Figure 1A). In rivers with dams or weirs downstream, water was sampled from both the upper and lower sides of the dams or weirs. Most of the sites were in the freshwater area, but some were influenced by the tide. The sampling was conducted in the warm season when eels are active.
The researchers who collected the water samples wore new disposable rubber gloves at each sampling site. Samples for the eDNA analysis were directly collected by suctioning the surface water 2-10 times (50 mL at a time), for which a disposable plastic pump (TERUMO syringe pump, Terumo Corp., Tokyo, Japan) was used. A plastic bucket fastened to a rope was used for water collection at some sites that could not be accessed. The bucket and tip of the rope tied to the bucket were decontaminated by spraying with sodium hypochlorite. Then, the bucket and rope tip were washed twice with environmental water. Samples measuring 50-600 mL of water were collected depending on the turbidity and immediately filtered at the sampling site using a syringe and an enclosed filter cartridge with a nominal pore size of 0.45 µm (Sterivex, Merck Millipore, Burlington, MA). When we could not filter immediately after sampling, benzalkonium chloride was added to the sampled water to a final concentration of 0.1% to suppress the degradation of DNA . The water samples were transported and filtered in the laboratory within 24 h. A 500 mL sample of pure water was filtered to monitor contamination at the filtration and subsequent extraction steps. Filters were preserved at −30 • C until DNA extraction.
Water temperature and salinity were measured (Horiba Kyoto, Japan, or YSI Inc., Ohio, United States) before the sample were collected, and the area surrounding the sampling points were surveyed for information on the environment including the presence/absence of a weir on the upper or lower side of the sampling point, bottom material, vegetation, and revetment.

eDNA Extractions and qPCR Analyses
The eDNA extractions and qPCR analyses followed the Environmental DNA Sampling and Experiment Manual (Minamoto et al., 2020). Before the commencement of the experiments, all instruments and equipment were cleaned with DNA-Off (Takara Bio Inc., Shiga, Japan) to avoid DNA crosscontamination. All extractions and quantitative polymerase chain reaction (qPCR) setups were performed in a designated eDNA laboratory in a room free of PCR products (i.e., no PCR machines and no prior PCR amplification occurring in the room). eDNA was extracted from the filters using the modified Qiagen DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany) extraction method (Miya et al., 2016). Briefly, the RNAlater solution was first removed by passing the Milli-Q water through the filter membrane. Then, 200 mL lysis buffer (AL) and 20 µL proteinase K were added to the filter, and the filter was incubated with rotation (20 rpm) (Roller 6 digital, IKA, Staufen im Breisgau, Germany) at 56 • C for 20 min. Subsequent extraction steps followed the standard Qiagen DNeasy extraction protocol. DNA was diluted using 100 µL elution buffer (AE), and all extracted DNA samples were preserved at −20 • C before qPCR experiments. Buffer AL, AE, and proteinase K were provided by the extraction kit.
Real-time quantitative PCR (qPCR) was conducted with a total reaction volume of 15 µL, which included 900 nM of each primer, 125 nM of a probe, 0.075 µL of AmpErase TM Uracil N-Glycosylase (Thermo Fisher Scientific, Waltham, MA, United States), and 2 µL of DNA template in 1x Environmental Master Mix 2.0 (Life Technologies, CA, United States) in a StepOne TM real-time PCR system (Life Technologies). Species-specific primers and a probe designed by Kasai et al. (2020) were used as follows: forward primer Aja-Dlp-F2 5 -TACATTTAATGGAAAACAAGCATAAGCC-3 , reverse primer Aja-Dlp-R3 5 -CGTTAACATTACTCTGTCAACTTACCTG-3 , and a probe Aja-Dlp-Pr 5 -FAM-ACCCATAAACTGATAAATAG-NFQMGB-3 . The primer set was specifically designed to discriminate A. japonica among the Anguilla species potentially inhabiting the inland and coastal areas of Japan because some foreign species such as the European eel (Anguilla anguilla) was imported and cultured in Japan (Kishikawa, 1997) and have been found in natural waters as they escaped from cultured ponds (Zhang et al., 1999). In Kasai et al. (2020), all of the other species and subspecies in the genus Anguilla were taken into account to design the assay, and in vitro test were also conducted to verify the specificity using mitochondrial DNA samples of A. anguilla, Anguilla bicolor bicolor, Anguilla marmorata, and Anguilla rostrata. Thermal cycle settings for PCR were as follows: 2 min at 50 • C, 10 min at 95 • C, and 55 cycles of 15 s at 95 • C, and 1 min at 60 • C. PCR was performed in triplicate for each eDNA sample. The limit for quantification was defined as three copies per reaction, and any one of the detections out of the three replicates was considered a true detection (Takeuchi et al., 2019;Kasai et al., 2020). For the PCR negative control, ultrapure water was used instead of the DNA template. Each PCR plate included a fourfold dilution series of commercially synthesized cloned DNA, including the target sequence of A. japonica (30,000 copies/reaction to 30 copies/reaction) and PCR negative control including 2 µL of ultrapure water in place of the template DNA, in triplicate.

Numerical Model of Larval Transport
Larvae of A. japonica are transported to Japanese coastal waters by the Kuroshio Current ( Figure 1B, Shinoda et al., 2011). Therefore, the distribution of larvae in terrestrial waters would depend on how they are transported by ocean currents and reach each river estuary. The swimming speed of larvae is slow compared to the velocity of the surrounding currents; thus, they are assumed to be particles floating in the water. If the details of the currents are recognized around Japan, it is possible to simulate how many larvae will reach each river.
Recently, the Japan Meteorological Agency developed a finescale ocean model (JPN model) covering the coastal seas of East Asia including Japan (117 • E-160 • E, 20 • N-50 • N, Sakamoto et al., 2019). The model can realistically reproduce the distributions and seasonal variations of water temperature, salinity, and sea level. The horizontal resolution of the JPN model is 1/33 • in the zonal direction and 1/50 • in the meridional direction, corresponding to approximately 2 km, which is sufficient for the reproduction of currents in the coastal regions. Small-scale features in the coastal seas, such as fronts and mesoscale eddies, were also well simulated with high resolution. In this study, therefore, larval transport was simulated using the daily average flow reproduced by the JPN model to investigate the recruitment success of Japanese eel.
As an initial condition, 80,000 particles that were regarded as larvae were released in the east of Taiwan because all the larvae transported to the waters around Japan originate from the Kuroshio Current (Hsiung and Kimura, 2019, Figure 1B). It is known that A. japonica larvae show diel vertical migration, remaining in the upper surface waters at night and diving to deeper waters during the day (Otake et al., 1998). Therefore, the particles were set to have diurnal vertical movement from the surface to a depth of 100 m (deeper in the daytime and shallower at nighttime) in the model. The particles were transported by the currents reproduced by the JPN model, and their positions were calculated by the Runge-Kutta fourth-order method. The simulation was conducted from 1 December to 30 April, based on the information of larval catch in the Japanese rivers (Yoneta et al., 2019). The recruitment success was evaluated by the cumulative number of particles reaching the shore. Particles were considered to have reached the shore if they were in the grids within a 2 km (horizontal resolution of the JPN model) distance from the shore. The simulation was repeated nine times using the velocity fields from 2008 to 2017. The results of the simulation were compared with those from the eDNA detection.

Collection of Geographical, Water Quality, and Demographic Data
To investigate the relationship between the river environment and eDNA concentration, which is supposed to represent the fish biomass, we used water quality data from the Water Quality Survey of Public Water Areas dataset of the Ministry of the Environment of Japan. The time-averaged data in 2016 were collected from water quality monitoring sites throughout Japan for dissolved oxygen (DO), chemical oxygen demand (COD), suspended solids (SS), total nitrogen (TN), minimum pH, and maximum pH. Basin mesh data (ver. 2.1) were collected from the National Land Numerical Information database of the Ministry of Land, Infrastructure, Transport, and Tourism of Japan. Water quality data within 1 km of the eDNA sampling site were selected and used for statistical processing. The function "it falls inside" in ArcMap 10.5.1 (ESRI, Redlands, CA, United States) was used to link 123 points.
The distance of each eDNA sampling point from the river mouth and the width of the river at the sampling point were measured using Google Earth Pro (ver. 7.3.2.5776, Map data: Google, Digital Globe). Data on the number of farmed eel that were released were taken from the Census of Fisheries (2003, 2008, and 2013).

Statistical Analyses
The Tukey HSD test was used to compare the differences in eDNA concentrations among the regions.
Generalized linear models (GLMs) with a Gaussian distribution were developed to identify surrounding environmental factors influencing eDNA concentration. Sites with positive DNA concentrations were employed in the model to assess the effect of the environment on the survival and growth of larvae in each river after they were transported. Previous comparisons between visual observation of fish and eDNA concentration in rivers suggest that eDNA quantification can detect the order of magnitude variations in fish abundance (Doi et al., 2017;Maruyama et al., 2018). Therefore, the logarithmically converted eDNA concentration (>0) was used as the dependent variable. Prior to modeling, the collinearity of all explanatory variables (distance from the river mouth, river width, water temperature, TN, DO, COD, SS, minimum pH, and maximum pH) was evaluated using correlation values and variance inflation factors (VIFs). Since the explanatory variables used in the model need to be independent of each other, VIFs were used to examine the multicollinearity among the explanatory variables and to remove variables that are highly correlated. As all values of the VIFs were below 5, non-collinearity was assumed in this study (Zuur et al., 2009; Table 1).
The best model was selected using Akaike's information criteria (AIC). The best model for any candidate set applied to a given dataset was that with the lowest AIC value, and models with AIC < 2 were assumed to be reasonable alternatives to the best model, and thus were retained (Burnham and Anderson, 2010). All statistical analyses, GLM analyses, and graphics employed R, version 4.0.2 (R Core Team, 2020), and the lme4 and MuMIn packages (Bates et al., 2015;Bartoń, 2019).

Environmental DNA Concentration
Successful amplification of eDNA was quantified at 181 sampling sites, accounting for 49.6% of the total sites that were sampled  ( Figure 2). In the PCR analysis, no positives were detected in the negative controls. The eDNA of eels was detected in some sites in the NSI and WSJ, although the concentrations were not high (Figure 3). In WK, PMJ, and SIS, the ratio of detected sites was over 80%, and the eDNA concentrations were high, although the variations were large with no detection in some rivers. The three largest values (44,356.6 copies/L from the Asaragi River, 41,177.0 copies/L, and 38,407.1 copies/L from the Aku River) recorded across all sampling sites were detected in the PMJ. The southern part of PT showed large values, while the detected values were small or zero in the northern part. The eDNA was rarely detected in NSJ and HD. In these regions, the proportion of sites where eDNA was detected was small, and even when detected, their concentrations were very low. The average eDNA concentration was significantly different among regions (p < 0.001, Table 2). The concentrations in HD and NSJ were significantly lower, while those in SIS and PMJ were significantly higher than in the other areas. Figure 4 shows the results of the simulation of larval transport; the number of particles that reached each coastal area by April 30. The simulation results are generally well consistent with the eDNA results (Figure 2). A large number of particles reached the coasts in NSI, PMJ, WK, and WSJ every year, while low in SIS. Particles reached the coastal areas of PT in half of the years we examined. Even in those cases, they only reached the south of PT and were rare in the north. Further, particles were rarely transported beyond the Noto Peninsula, and they never reached HD. Only a few particles reached the east of the Noto FIGURE 3 | Concentrations of environmental DNA in each region. In the box plots, the midline represents the median, the upper and lower limits of the box represent the third and first quartiles, respectively, the whiskers represent 1.5 times the interquartile range from the top (bottom) of the box to the furthest datum within that distance, and the circles represent the individual data beyond that distance. The numbers given above the graph are the number of sites where eDNA was detected (DS), was not detected (NDS), and the ratio of detected sites to the total number of sites that were sampled in each region. The regions include Nansei-Shoto Islands ( The number of particles reaching an area was compared with whether eDNA was detected. GLM analysis was performed using a binomial distribution with eDNA detection as the dependent variable and the average number of particles reaching the area as the explanatory variable ( Figure 5). The binomial test showed that eDNA was detected in the rivers where the number of particles reaching the area was high (p < 0.01).

Important Environmental Factor for Eel Survival and Growth
To determine the river environments that are most suitable for eels, we modeled the detected eDNA concentrations using environmental factors as explanatory variables. Among the GLMs used to examine eDNA concentration, seven models with AIC < 2 were retained ( Table 3). Four common explanatory variables (distance from the river mouth, TN, DO, and the minimum pH) were selected in these models. The relationship of the eDNA concentration to the distance from the river mouth was negative, while that with TN, DO, and the minimum pH was positive. Only TN concentration was significantly correlated with eDNA concentration.

Effect of Released Eels on eDNA Detection
Young eels are released by prefectural Inland Fisheries Cooperatives throughout the country, other than in HD, to allow their growth and to increase their population sizes in the natural freshwater conditions. This indicates that individuals that were grown in eel farms and artificially released and those that had migrated to the natural environment as glass eels were mixed in the natural environment. Because eels form a single spawning population throughout their distribution area, it is impossible to genetically discriminate between naturally recruited and released individuals. Figure 6 shows the number of farmed eel individuals released in each prefecture. Large numbers were not released in any of the regions, but were released evenly throughout the country except for HD. There was no correlation between the number of released eels and eDNA (Figure 2, p > 0.05). Thus, although some of the results of the eDNA distribution in this study may

DISCUSSION
This is the first study to investigate the distribution of Japanese eels in detail on a national scale, and this could not be done without using the eDNA method. The eDNA of Japanese eels was detected in an unexpectedly wide area (Figure 2), even though A. japonica is listed as an endangered species by both the Ministry of the Environment and the IUCN. The detection rates were as high as over 80% in PMJ, WK, and SIS (Figure 3). This result coincides with the distribution of larval reaching rates estimated by the numerical simulation (Figures 4, 5, p < 0.01). This indicates that the major distribution of A. japonica in Japanese rivers is determined by the maritime transport of larvae. Although a quantitative relationship between eDNA concentrations and biomass has not been assured, it can be assumed that the eDNA concentrations are reflective of the biomass (Maruyama et al., 2018;Itakura et al., 2019). Generally, our results revealed that eels in the river are distributed south of Noto Peninsula, in the Sea of Japan, and south of Miyako on the Pacific side. This is roughly consistent with the results compiled and reviewed by Yoneta et al. (2019), although our results show a much more detailed distribution. These areas are presumed to be continuously populated by naturally recruited individuals based on the results of the larval transport simulation. This is also consistent with recent reports on the harvests of glass eels, the collection of Japanese eels in river systems where eels were not released, and otolith stable isotope ratio-based origin determinations (Yoneta et al., 2019). The small amounts of eDNA detected north of PT (Figure 2) possibly reflect the artificially released individuals, because only one of the 19 individuals was a natural individual from a lake in the north of PT (Yoneta et al., 2019). The distribution of eDNA showed a clear difference between the Sea of Japan side and the Pacific side of the Tohoku region (Figure 2). eDNA was detected in the southern region of PT rivers, while only local and very little eDNA was detected in the NSJ. This difference reflects the transport process of A. japonica to Japanese coastal waters from the south. Although some previous studies suggest the existence of natural individuals in PT (Kume et al., 2020), the transport process of larvae to PT rivers is unknown. In other words, it has been unclear whether eel larvae are transported to the south of PT by the Tsushima Current from the Sea of Japan through the Tsugaru Strait or migrate northward from the Kuroshio and/or Kuroshio Extension (Figure 1B). To address this problem, our eDNA distribution (Figure 2) and numerical experiments (Figure 4) indicate that the latter process is plausible. Generally, the Kuroshio flows eastward off Cape Inubozaki, but a part of it often moves northward as a branch along the Pacific coast of Japan ( Figure 1B). It is also well known that the Kuroshio meanders east of the Japanese Islands. As it meanders and is unstable, a part of the Kuroshio sometimes cuts off the mainstream and moves northward as a mesoscale eddy (Yasuda et al., 1992). Eel larvae would use the branch or warm core rings to migrate northward and reach the south of the PT coast.
From the perspective of the larval transport in the ocean, the occurrence of these eels should have been more in the NSI, but this was not reflected in our eDNA survey. There are two possible reasons for this discrepancy. The first reason is the shortness of rivers in the NSI. There are no large rivers, and the river discharges are small in small islands. Glass eels approach estuaries and rivers by smelling the odor originating from rivers (Fukuda et al., 2019). Therefore, around NSI, where there are no large rivers, eel larvae cannot smell the river and cannot approach the mouth of the river. The second reason is habitat segregation. Kano et al. (2017) surveyed the occurrence of A. japonica and the Indo-Pacific eel (A. marmorata) in many streams in Kyushu and NSI. They found a considerably higher abundance of A. marmorata than A. japonica in NSI, while no A. Marmorata was found in Kyushu. A. marmorata is widely distributed in tropical and subtropical terrestrial waters (Watanabe et al., 2008). A. japonica may be outcompeted by A. marmorata, which is larger than A. japonica, and their abundance is hence likely to be low in NSI. For these reasons, Japanese eels may not be well adapted to living on small islands where a tropical anguillid eel resides.
Based on our simulation, only a few particles reached the SIS; however, this is not the true picture as shown by the eDNA detected. This divergence is likely due to the characteristics of the model. As shown in Figure 1B, the two entrances of SIS, the Bungo Channel and the Kii Channel are very narrow with complicated topographies, making the friction and viscosity large in the model. In addition, it is difficult to reproduce the sufficient water exchange through these channels, because the minimum scale of the phenomena that can be reproduced by the model is about 5-10 times the grid size of the model (Lévy et al., 2012). Therefore, it is difficult for particles to pass through the channels and enter the SIS. Another improvement is the preference for low salinity areas. Eel larvae detect odor and migrate to estuaries. Including the effect of low salinity preference into the model may help to reproduce a more realistic distribution of eDNA.
The distance from the river mouth, TN, DO, and the minimum pH were selected in the GLMs. Considering the fact that the larvae migrate from the sea to the river, it is a natural consequence of the negative correlation between the distance from the river mouth and the eDNA concentration. Since oxygen is necessary for all animals to sustain life, it is not surprising that there was a positive correlation between DO and eDNA concentrations. Rather, it is worth noting that TN was a significantly important factor for determining the concentration of eDNA (Table 3). This indicates that the survival and growth of eels are better in nutritious rivers in Japan. Japanese rivers were highly eutrophicated in the 1970s and the 1980s, when there was a large growth in the Japanese economy. Red tides and consequent hypoxia often occurred, and the biodiversity in the rivers was reduced during this high economic growth period. Since the 1990s, however, sewage treatment technology has been developed, and the inputs of nutrients and organic matter from households and industries have been reduced to protect the environments of rivers, estuaries, and coastal areas. This policy has worked, and the river water condition has drastically recovered (Ye and Kameyama, 2020), and many species are returning to the Japanese rivers. The average TN value in this study was 1.2 mg/L. Considering that TN concentrations in formerly eutrophic rivers exceeded 5 mg/L in the 1980s (Ye and Kameyama, 2020), this value is very low. Therefore, the current high TN value does not necessarily mean that the river is highly eutrophicated. The high abundance of eels in TN-rich rivers is likely to be related to the abundance of food.
Unlike conventional methods that capture fish, the eDNA method requires less labor and a shorter time and reduces onsite survey costs. It enables surveys to be conducted over a wider area at multiple locations. In addition, unlike conventional surveys, the eDNA survey does not use fishing gears, and no special techniques are required to collect samples and only a small amount of water needs to be collected. Thus, these surveys can be conducted easily. This makes it possible to conduct unprecedented multisite surveys, similar to that undertaken in this study. Such a multisite survey in a short period would not be possible using conventional methods. Furthermore, this method does not kill or injure any organisms by capturing the target species and requires little access to their habitats. It is especially beneficial for endangered and rare species such as the Japanese eels.
One drawback of the eDNA method is that the demography of the population cannot be determined. The length and weight of the individuals of the target species also cannot be estimated. In the case of eels, we also cannot determine whether the individuals originated from farms or from the wild. Further, the eDNA also flows down from the upstream side of the river, and hence may not reflect the biomass at the site, but the sum of the transported and on-site DNA (Jerde et al., 2011). Therefore, if necessary, a combination of eDNA with conventional capture surveys might be helpful to obtain more detailed information that is needed for conservation planning and to take prompt and effective measures.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.