Distinct seasonal migration patterns of Japanese native and non‐native genotypes of common carp estimated by environmental DNA

Abstract Understanding behavioral differences between intraspecific genotypes of aquatic animals is challenging because we cannot directly observe the animals underwater or visually distinguish morphologically similar counterparts. Here, we tested a new monitoring tool that uses environmental DNA (eDNA), an assemblage of DNA in environmental water, to specifically detect Japanese native and introduced non‐native genotypes of common carp (Cyprinus carpio) in Lake Biwa, Japan, and estimated differences between the two genotypes in the use of inland habitats. We monitored the ratios of native and non‐native single nucleotide polymorphism alleles of a mitochondrial locus of common carp in a lagoon connected to Lake Biwa for 3 years using eDNA. We observed seasonal dynamics in the allele frequency showing that the native genotype frequency peaked every spring, suggesting that native individuals migrated to the lagoon for spawning and then returned to the main lake, whereas non‐native individuals tended to stay in the lagoon. The estimated migration patterns corresponded with the estimates of a previous study, which were based on commercial fish catch data. Our findings suggest that eDNA‐based monitoring can be useful tool for addressing intraspecific behavioral differences underwater.


| INTRODUCTION
Intraspecific behavioral differences driven by genetic variation are an important factor that affects ecological and evolutionary processes (Wolf & Weissing, 2012). However, in situ observation of such differences is generally very difficult in aquatic animals because of the inability to observe underwater animals and the difficulty in visually discriminating morphologically similar counterparts. A new development in monitoring animal distribution during the last decade can fill this methodological gap, as this method uses environmental DNA (eDNA) as a trace for animal presence. eDNA is defined as DNA assemblages existing in environmental media, such as water and soil. As genetic assessment was first applied to eDNA recovered from water to estimate the distribution of bullfrogs (Ficetola, Miaud, Pompanon, & Taberlet, 2008), eDNAbased monitoring of animal distribution has rapidly developed and gained general acceptance for use in detecting specific species (see reviews by Lodge et al., 2012;Rees, Maddison, Middleditch, Patmore, & Gough, 2014;Thomsen & Willerslev, 2015) and taxa (Miya et al., 2015;Thomsen et al., 2012;Valentini et al., 2016) in aquatic systems.
The eDNA method can precisely identify species from genetic information, which is especially useful when dealing with similar-looking organisms, such as intraspecific genotypes. We previously developed an eDNA method to quantify the relative biomass ratio of intraspecific genotypes of common carp (Cyprinus carpio) based on a single nucleotide polymorphism (SNP) (Uchii, Doi, & Minamoto, 2016). Here, we tested the potential use of this method to estimate the difference in the seasonal migration patterns of two genotypes of common carp.
As the Japanese archipelago broke apart from the Eurasian continent, a unique strain of common carp has evolved in Japan. Based on mitochondrial sequence data, the Japanese strain is estimated to have diverged from Eurasian ones ca. two million years ago (Mabuchi, Senou, Suzuki, & Nishida, 2005). However, domesticated Eurasian strains were introduced to Japan in 1905 at the latest (Maruyama, Fujii, Kijima, & Maeda, 1987), resulting in today's extensive invasion of non-native common carp throughout Japanese freshwater systems (Mabuchi, Senou, & Nishida, 2008). Lake Biwa, the largest and oldest lake in Japan, is an important habitat for the Japanese native strain and other endemic freshwater fishes of Japan. An early report suggested that the native strain inhabited the main lake and migrated to littoral and inland habitats in the spring for spawning, whereas the nonnative strains inhabited the littoral and inland habitats all year round (Hurukawa, 1966). This estimate was based on only 1 year of commercial catch data in which the two strains were visually discriminated; thus, we still do not have concrete evidence about the migration patterns of native and non-native strains.
In this study, we applied eDNA monitoring to estimate the seasonal use of inland habitats by the native and non-native genotypes.
The ratio of the two genotypes in a lagoon connected to Lake Biwa was monitored for 3 years using a previously developed eDNA assay that quantifies the ratio of native and non-native DNA alleles based on a SNP (Uchii et al., 2016). We observed reproducible seasonal dynamics in the ratio of two genotypes, which suggested the differences in the migration patterns between the two genotypes.

| Water sampling and filtration for eDNA collection
Water sampling was performed in Ibanaiko Lagoon (Fig. 1), which is connected to Lake Biwa, Japan, 19 times between April 2013 and April 2016. To cover the wide range of the lagoon system, we set four monitoring sites on the shore, three of which were located from upstream to downstream in the lagoon (sites a-c; Fig. 1), and one of which was located in a channel flowing from the lagoon to the main lake (site d; Fig. 1). We collected 2 L of surface water at the four shore sites during every sampling occasion. In total, 76 samples (4 sites × 19 sampling times; Table 1) were collected on the shore. We also collected 2 L of surface water from three offshore sites (sites e-g; Fig. 1 with detergent and distilled water in the laboratory, were rinsed with environmental water collected on site at least three times before sampling. Water samples were immediately placed in a cold box and transported to the laboratory within 5 hr. Each water sample was filtered onto a Whatman ® GF/F filter (47-mm diameter, 0.7μm particle retention size; cat no. 1825-04, GE Healthcare Japan, Tokyo, Japan) to capture eDNA. The filtered volume per filter ranged from 330 ml to 1 L but was 500 ml in most samples (see Table S1). The filtration apparatus was decontaminated with 5% bleach for at least 5 min and thoroughly washed with tap water and distilled water between samples. At the end of each sampling day, 1 L of distilled water was subjected to filtration in the same manner as the samples to control cross-contamination during filtration and subsequent DNA extraction procedures. The filters were stored at −20°C until DNA extraction.

| eDNA extraction from GF/F filters
DNA was extracted from each GF/F filter using the DNeasy ® Blood & Tissue Kit (Qiagen, Hilden, Germany). Each filter was soaked with 400 μl of Buffer AL and 40 μl of Proteinase K in an internal container of a Salivette ® tube (Cat no. SAR-511534; Sarstedt, Nümbrecht, Germany) and incubated at 56°C for 30 min. After centrifugation at 5,000 g for 5 min, 220 μl of TE buffer (pH 8.0) was added to each filter in the tube, and the tubes were centrifuged again under the same conditions. Buffer AL (200 μl) and 100% EtOH (600 μl) were then added to each filtrate and mixed by pipetting. The mixture (~600 μl) was applied to a DNeasy Mini Spin Column and centrifuged at 6,000 g F I G U R E 1 Sampling locations in Ibanaiko Lagoon, which is connected to Lake Biwa. Shore sites are labeled a-d and offshore sites are labeled e-g Offshore sites 0 0 4 3 for 1 min. This step was repeated until the mixture was completely processed. We followed the manufacturer's instructions for further steps. Finally, DNA was eluted from the column with 100 μl of Buffer AE, with some exceptions (150 μl; Table S1).

| Quantification of ratios of native and nonnative alleles in eDNA
Cycling probe technology is a highly sensitive method that uses DNA-RNA-DNA chimeric probes to detect SNPs; it involves using ribonuclease H to cut the RNA portion of the chimeric probe that hybridizes to the target DNA, thereby enabling sequence-specific detection (Bekkaoui, Poisson, Crosby, Cloney, & Duck, 1996;Yatabe et al., 2006). We quantified the relative abundances of the DNA al- The last base (A) at the 3′ end of the original reverse primer (Uchii et al., 2016) was removed to increase mismatches to related species within the last five bases at the 3′ end (Appendix 1). The standard mixtures of native and non-native DNA (50-100 copies/reaction) in the ratios of 15:1, 7:1, 3:1, 1:1, 1:3, 1:7, and 1:15 (15, 7, 3, 1, 0.33, 0.14, and 0.07 in terms of native/non-native DNA concentration ratios, respectively) were prepared, and at least five of them were included in triplicate or duplicate in every run (see details in Table S2). Using standard curves created between ΔC T (calculated as C T(non-native DNA) -C T(native DNA) ) and the native/non-native DNA concentration ratios of the standard mixtures, the native/non-native DNA concentration ratios of eDNA samples were quantified when the amplification signals were detected at C T < 40 for both native and non-native probes in at least two PCR replicates of a sample. A replicate in which either probe showed no signal was omitted from the calculation. At least two PCR blanks were included in every run to control for DNA contamination in the PCR. See Appendices 2 and 3 for the PCR inhibition test.

| Quantification of common carp cytochrome b DNA in eDNA
When estimating the native and non-native allele frequencies in the lagoon, we quantified the amount of common carp cytochrome b DNA (CytB), which was previously demonstrated to be linearly correlated with carp biomass in aquarium experiments (Takahara, Minamoto, Yamanaka, Doi, & Kawabata, 2012), to adjust for variations among sites in the amount of DNA. The copy numbers of CytB were quantified according to Takahara et al. (2012). Briefly, real-time PCR was performed in triplicate in a reaction volume of 20 μl containing 1× TaqMan Table S3). The CytB concentrations of eDNA samples were quantified when the amplification signals were detected in at least two PCR replicates of a sample. The concentration of CytB at each site (copies/L) was calculated from the volumes of filtered water (330-1,000 ml) and eluted DNA (100-150 μl).

| Estimation of native and non-native genotype frequencies using eDNA
The ratio of native DNA was calculated for each eDNA sample as fol- The time-series data of the genotype frequencies for 3 years at each of the four shore sites were checked using a unit root test to discard the possibility of randomness. We performed the Phillips-Perron unit root test for the null hypothesis that the time series has a unit root against a stationary alternative. To check the possibility that the genotype frequencies were linearly decreased/increased over time, the genotype frequencies of the four shore sites were also analyzed by a binomial generalized linear model (GLM) with logit link function, with sampling year and site as explanatory factors. All statistical analyses were performed using R ver. 3.3.1 (R Core Team 2016),

| Contamination during filtration, DNA extraction, and real-time PCR
No amplification signals were detected in the eDNA extracts from the control GF/F filters used to filter 1 L of distilled water in the real-time PCRs quantifying the native/non-native DNA ratios and CytB. No amplification signals were detected in the PCR blanks in any real-time PCR assays.

| Dynamics of the frequencies of native and nonnative genotypes
The native/non-native DNA concentration ratios, which were quantified based on the SNP in the D-loop region, were successfully estimated for 85 of 97 samples (see Table S2 for C T values and standard curves). The 12 eDNA samples that failed to meet estimates had low copy numbers of common carp mtDNA (i.e., <1 copy/μl of CytB). or interactions between these two factors (Chi-square = 1.12, df = 9, p > .9), discarding the possibility that the genotype frequencies were linearly decreased/increased over time. The nonweighted frequencies of the native genotype estimated by a simple averaging of the ratios of native DNA at the shore sites showed repetitive seasonal dynamics in which the frequencies were highest in the spring when spawning occurred and lower in other seasons (Fig. 2a, plotted as open squares). The nonweighted frequencies of the native genotype estimated from the shore and offshore sites combined (Fig. 2a, open circles) showed similar values to those estimated from the shore sites alone (open squares). The CytB-weighted frequencies of the native genotype ( Fig. 2b) showed slightly different dynamics compared with the nonweighted frequencies. However, the overall trend was similar, as the CytB-weighted frequencies of the native genotype were also highest in the spring every year. The CytB-weighted frequencies estimated from the shore sites alone (Fig. 2b, solid squares) showed similar values to those estimated from the shore and offshore sites combined (solid circles).

| DISCUSSION
We estimated the differences in the seasonal use of inland habitats between Japanese native and introduced non-native genotypes of common carp in Lake Biwa by monitoring SNP allele frequencies in eDNA. Samples collected over 3 years, including four springs, showed repetitive seasonal dynamics in the frequency of the native genotype for both the nonweighted (Fig. 2a, open squares) and CytB-weighted These frequencies were calculated when the ratios of native DNA were quantified at more than two sites on each sampling day. Vertical bars represent standard errors. The ratios of native DNA at each site (a-g) are also plotted. (b) CytB-weighted frequencies of the native genotype in the lagoon estimated from the shore sites (solid squares) and the shore and offshore sites combined (solid circles). These frequencies were calculated when both the ratios of native DNA and the concentrations of CytB were quantified at more than two sites on each sampling day suggesting that many individuals with the native genotype left the lagoon and moved, presumably to the main lake, whereas non-native individuals tended to remain in the lagoon. Such a difference in migration patterns between the native and non-native genotypes corresponds well to the previous estimate based on commercial fish catch data (Hurukawa, 1966). The latter showed that the catch amount of the native strain increased in the spring and decreased after the summer in littoral and inland habitats and increased in the winter in pelagic habitats. However, the catch amount of the non-native strains did not fluctuate substantially in the littoral and inland habitats and was very low in the pelagic habitats throughout the year (Hurukawa, 1966). Furthermore, the seasonal changes in the native genotype frequency observed in the present study were roughly consistent with the changes detected in common carp individuals that were caught in Ibanaiko lagoon in a previous study (Uchii, Okuda, Minamoto, & Kawabata, 2013; see Appendix 4). These observations suggest the potential use of eDNA for detecting behavioral differences between closely related species.
When we deal with morphologically similar organisms, such as intraspecific genotypes and subspecies, there is often difficulty in visual discrimination. The eDNA method can readily discriminate genotypes because it uses genetic information. Furthermore, we can obtain long-term data of genotype frequencies with substantially less effort than is needed for traditional catch sampling or biotelemetry (i.e., genetic analyses of large numbers of individual fish after catch or long-term tracking of large numbers of individuals after tagging them, respectively). As research generally requires long-term observations to detect seasonal patterns of animal migration, the time-and costeffectiveness of the eDNA method would be a great advantage.
On small spatial or short temporal scales, there is difficulty in linking eDNA information with animal presence and biomass due to eDNA dispersal, eDNA degradation, and animal movement (Barnes & Turner, 2016;Barnes et al., 2014). Thus, eDNA monitoring would be more suitable for tracking animal movement at larger spatial and temporal scales, such as seasonal migration, because the dispersal range of eDNA and the stability of DNA molecules are limited in the environment (e.g., Barnes et al., 2014;Jane et al., 2015;Maruyama, Nakamura, Yamanaka, Kondoh, & Minamoto, 2014;Strickler et al., 2015). A recent study demonstrated that the existence of target DNA was correlated with the migration range of anadromous fish (Yamanaka & Minamoto, 2016). Another study demonstrated a significant positive correlation between DNA concentration and the number of bighead carp detected by telemetry in a spawning habitat in which the increase in DNA was presumably attributable to the individuals that migrated for spawning (Erickson et al., 2016). These studies, together with the present study, suggest the great potential of the eDNA method for monitoring seasonal animal movement.
Most eDNA studies conducted thus far have used mtDNA markers to increase the detection probability because the copy number of mtDNA is much greater than that of nuclear DNA. We also targeted mtDNA for the same reason. However, mtDNA markers have a disadvantage in that they do not distinguish hybrids due to the maternal inheritance of mitochondria. As intraspecific genotypes can hybridize, we cannot exclude the possibility that the native genotype frequency based on mtDNA haplotypes would be over-or underestimated if cytonuclear disequilibrium was to exist. The use of nuclear DNA markers would solve this problem. Although several nuclear DNA markers that distinguish the Japanese native and non-native genotypes of common carp have been reported (Mabuchi, Song, Takeshima, & Nishida, 2012), these markers are single-copy and thus very difficult to detect in eDNA. However, because these markers can detect intraspecific genetic structure, the development of an effective eDNA concentration method to detect single-copy genes would greatly expand the scope of eDNA studies. Alternatively, the development of markers in multiple-copy nuclear genes would offer great potential for evaluating intraspecific genetic structure. Nuclear ribosomal DNA (rDNA) would be a strong candidate, as a recent study demonstrated that a marker in the internal transcribed spacer region of rDNA has greater sensitivity than a mtDNA marker in eDNA detection . As many eDNA studies have repeatedly noted, the greatest advantage of eDNA monitoring is the simplicity of its sampling and analysis pro-

O)
A P P E N D I X 2 Shift of threshold cycle (C T ) values in field eDNA samples spiked with the standard DNA mixture of native and non-native DNA at a ratio of 1:1 (100 copies each/reaction). The copy number of the spiked DNA was ca. >10 times larger than that of target the DNA originally contained in the filed samples. Real-time PCR was performed in duplicate. Delayed C T values were observed in some samples when 5 μl of eDNA was used as template, whereas little or no inhibition effect was detected with 2 μl of template. Non-native A P P E N D I X 4 Frequency of the native genotype by month for common carp (≥350 mm in standard length and appeared to be adult) captured in Ibanaiko Lagoon in 2006, and 2009in Uchii et al. (2013. Because common carp are no longer commercially fished in this area, all carp were caught as bycatch. Note that the native strain of common carp is more susceptible to the lethal pathogen Cyprinid herpesvirus 3, which was introduced to Lake Biwa in 2003, than are the non-native strains (Ito, Kurita, & Yuasa, 2014); thus, the proportion of the native strain in Lake Biwa is likely to gradually decline with time (Uchii et al., 2013)