Introduction

How genetic structure is shaped across a landscape is an essential theme in evolutionary biology, and this information provides an empirical basis for conservation biology by elucidating habitat connectivity or predicting the effects of landscape change (Hohenlohe et al. 2021). A pattern of ‘isolation by distance’ (IBD; Wright 1943), in which genetic differentiation increases with geographic distance, is a very common pattern to explain population structure (Meirmans 2012). In addition to space, another important component of a landscape is the environment (Nosil et al. 2005; Thorpe et al. 2008). In the field of landscape genetics, many studies have investigated the interactions between genetic differentiation and landscape variables, and a couple of patterns have described the relationship between genetic structure and environment. When gene flow between populations that inhabit different environments is limited due to local adaptation or other factors, genetic differentiation increases with environmental differences and a pattern of ‘isolation by environment’ (IBE; Wang and Summers 2010) is generated. Another pattern that involves genetic structure and the environment is ‘counter-gradient gene flow’ (Sexton et al. 2014). In this scenario, gene exchange between different environments is high and forms a source-sink pattern in a landscape. Because these interactions between genetic differentiation and landscape variables are critically important for addressing classical evolutionary questions related to ecological speciation or conservation, the environment is a factor that cannot be ignored when describing the observed population structure (Orsini et al. 2013; Wang and Bradburd 2014).

However, it is very difficult to disentangle the relative strength of space (distance) and environment in observed patterns of spatial genetic differentiation because environmental differences are usually highly correlated with geographic distance (Sexton et al. 2014; Wang and Bradburd 2014). To test spatial relationships, the Mantel test and partial Mantel test are commonly used. However, the Mantel test family is strongly criticized for its inflated Type I error rate, especially under conditions where the measured variables are spatially correlated (Guillot and Rousset 2013; Harmon and Glor 2010; Meirmans 2012; Raufaste and Francois 2001; Rousset 2002). Although there are several proposed alternative methods to the Mantel test such as multiple regression on distance matrices (MRM; Lichstein 2007), Procrustes analysis, or redundancy analysis (RDA), and some of them perform much better than the partial Mantel test, previous studies have suggested that all methods fail to correctly model the relative importance of space and environment under certain patterns, which cannot be known a priori (Diniz-Filho et al. 2013; Gilbert and Bennett 2010). It is now recognized that no method has the ability to address the spatial-environmental correlation and perform hypothesis testing simultaneously (Zeller et al. 2016). These studies suggest the difficulty of separating spatial and environmental effects on genetic differences. However, once a sampling design that decouples distance and environment is developed, we will be able to clearly compare the relative strength of space and environment independent of these controversies. In the field of community ecology, Gilbert and Lechowicz (2004) used a sampling design that removed spatial autocorrelation of the environment sampled, and they analyzed the species’ spatial and environmental correlations to reject the neutral theory of biodiversity. However, subsequent attempts to partition spatial and environmental components have focused exclusively on statistical and analytical methodology (Gilbert and Bennett 2010). Especially in landscape genetics, sampling design has not received much attention (Meirmans 2015).

Here, we focused on stream ecosystems, which often exhibit heterogeneous environments (Heino et al. 2013; Uno 2016). In stream ecosystems, water temperature is the strongest factor affecting river organisms as well as flow regime (Olden and Naiman 2010), and it is increasingly being an important variable under ongoing climate change. Water temperature is determined by air temperature, riparian forests, and groundwater, among others (Caissie 2006; Nakamura and Yamada 2005). In particular, the effect of groundwater discharge on water temperature has been indicated to be very strong (Arscott et al. 2001). However, the influence of groundwater discharge on the upper sections of streams has received much less attention (Brown et al. 2007), and many ecological studies have used air temperature instead of water temperature (e.g., Almodóvar et al. 2012; Middaugh et al. 2018). Groundwater shows a very stable water temperature and flow regime; hence, groundwater discharge greatly affects the stream environment (Poff et al. 2010). Given water temperature in summer that is critical for cold-water organisms, streams with high groundwater display lower values and these streams could be an important habitat for these species. Since the amount of cold groundwater discharge vary locally within a watershed, groundwater creates spatial heterogeneity of environments within the watershed. Thus, to reveal the ecological role of water temperature, knowledge of the spatial heterogeneity caused by groundwater is thought to be useful.

As a practical study, we investigated the genetic structure of the cold-water fluvial sculpin Cottus nozawae in the upstream section of the Sorachi River, Hokkaido, Japan. Cold-water fish are very vulnerable to climate change, and it is worth studying the relationships between population structure and water temperature. Together with C. nozawae, Salvelinus malma, Salvelinus leucomaenis, Barbatula oreas, and Parahucho perryi inhabit this watershed. However, studies on S. malma and S. leucomaenis did not display strong genetic structure within this watershed due to their active migration (Koizumi 2011; Nakajima et al. 2020). Due to the lower mobility of adult C. nozawae than these species (Goto 1998; Okumura and Goto 1996), we predicted that C. nozawae should display clear genetic differentiation within the watershed. In addition, the summer water temperature has been shown to be the dominant factor determining the distribution and survival of this species (Yagami and Goto 2000), and a relationship in which low-water temperature in summer increases the density of this species has also been confirmed in our study area (Suzuki et al. 2021). Therefore, summer water temperature can be explicitly used as an important variable for C. nozawae. Hypothesis-driven studies that focus on a given variable have advantages over exploratory studies in designing sampling strategies (Richardson et al. 2016). Since ecologically important variables affecting survival and population density could cause local adaptation between different environments (Kawecki and Ebert 2004), we predicted that the differences in summer water temperature cause local adaptations and lead to the IBE pattern within the watershed.

The aims of this study are (i) to minimize the correlation of geographic distance and water temperature differences via a groundwater-focused sampling strategy, (ii) to investigate the local genetic structure of C. nozawae and its determinant factors, and (iii) to discuss the relationship between water temperature and the population structure of cold-water fish. Although recent studies have frequently used adaptative genetic markers to detect local adaptation, barriers to gene flow imposed by selection and local adaptation between populations can be detected with neutral markers (Sexton et al. 2014). Given that the use of adaptive markers is still challenging, the evaluation of current genetic differentiation and connectivity based on neutral genetic markers is still informative and can assist in allocating conservation units to preserve local genetic variation (Tsuda et al. 2015). To assess the genetic structure, we used putatively neutral genome-wide SNPs on the inter simple sequence repeat (ISSR) region and analyzed genetic differentiation and gene flow patterns.

Materials and methods

Study sites and sampling

This study was conducted in the upper section of the Sorachi River, Hokkaido, Japan (Fig. 1; Table 1). To conduct sampling across a heterogeneous landscape, we used an approach based on watershed geology, which is known to cause significant variations in water temperature through groundwater discharge (Caissie 2006; Nagasaka and Sugiyama 2010; Tague et al. 2007). In the Sorachi River, many small spring-fed inputs are found in the Quaternary volcanic region (Koizumi and Maekawa 2004; Watz et al. 2019). Thus, the sampling sites were selected so that tributaries with volcanic watersheds and nonvolcanic watersheds were spatially intermingled. We also chose sites interspersed throughout the watershed such that large spatial gaps among populations were eliminated as much as possible. Watershed geology was assessed based on the Seamless Digital Geological Map of Japan V2 from the National Institute of Advanced Industrial Science and Technology. There is a large dam (Kanayama Dam) in the watershed. While the area upstream of Kanayama Dam represents one of the most continuous river habitats in Japan and there are no barriers between populations, the downstream region is slightly more influenced by anthropogenic impacts on streams. For example, one population (Pop13) is located upstream of a check dam, and another population (Pop19) is located in a tributary with a downstream portion that runs through an agricultural landscape.

Fig. 1: Sampling localities and population structure.
figure 1

A Study area. The fill colors of the sampling site depend on the maximum water temperature. Site labels correspond to population IDs in the text and Table 1. B Genetic structure inferred by STRUCTURE. The number indicates the population IDs.

Table 1 Details of the sampling sites and genetic diversity of each population.

From 20 sites in the Sorachi River, 531 individuals of C. nozawae were caught for genetic analysis. The fish were caught by electrofishing (model 12-B Backpack Electrofisher; Smith-Root Inc.) or with hand nets. Small pieces of fin tissue were clipped, placed in 99.5% ethanol, and stored at −20 °C in the laboratory until DNA extraction. Total DNA was extracted using the QIAGEN DNeasy Blood and Tissue Kit (QIAGEN Inc.) with a combination of Genomic DNA Extraction Column (FAVORGEN Inc.).

Water temperature

At each sampling site, the summer water temperature was recorded hourly using data loggers (Onset Computer Corp. or Gemini Data Loggers Ltd., depending on site). Water temperatures were measured during the summer of 2020, and also during the summer of 2019 at some locations. To align the data sampling periods, data from 16 July 2020 to 31 August 2020 were extracted, and the maximum water temperature in this period was identified. At one site (Pop10) where we failed to record water temperatures in some mid-summer periods, the maximum temperature in 2019 was used as an alternative. We used the maximum water temperature which significantly affects the populations of cold-water fish as a representative of summer water temperature (Yagami and Goto 2000). The correlation of water temperature differences with waterway geographic distance (hereafter referred to as geographic distance) was examined with the Mantel test (9999 permutations) using the package ECODIST 2.0.3 (Goslee and Urban 2007) in R 3.6.0 (R Core Team 2019). Then, a Mantel correlogram showing the spatial correlation of water temperature differences across multiple ranges of geographic distance was constructed. The Mantel correlogram in ten equal distance classes was assessed with 9999 permutations using the package VEGAN 2.5.6 (Oksanen et al. 2019) in R.

MIG-seq library preparation and sequencing

To obtain neutral genome-wide SNP data, we used multiplexed ISSR genotyping by sequencing (MIG-seq; Suyama and Matsuki 2015), a technique in which loci between two microsatellite regions are amplified by PCR and next-generation sequencing. MIG-seq is one of the reduced representation sequencing methods along with restriction site-associated DNA sequencing (RAD-seq), but the number of available informative loci (MIG-loci) detected by MIG-seq is less than that of RAD-seq, and most of the MIG-loci are putatively neutral. Although MIG-seq is not suitable for outlier analysis or gene identification, acquired loci are sufficient for population genetic analysis. A MIG-seq library was prepared following the protocol outlined in Suyama and Matsuki (2015), except for the minor modifications outlined below. The first PCR was conducted using eight ISSR primer sets with tail sequences at an annealing temperature of 38 °C. The second PCR was conducted using primer pairs including tail sequences, adapter sequences for Illumina sequencing, and five-base (forward) and nine-base (reverse) barcode sequences to identify each individual sample. The conditions for the second PCR were as follows: 12 cycles of denaturation at 98 °C for 10 s, annealing at 54 °C for 15 s, and extension at 68 °C for 1 min. The second PCR products of all 531 samples were mixed, and then fragments in the size range of 400–800 bp in the purified library were isolated. After library quantification, the products were sequenced on the Illumina MiSeq platform using MiSeq Reagent Kit v3 (150 cycles). Both ends of the fragments and the indexes were read by paired-end and index sequencing: 80, 80, 9, and 5 bases of sequences were determined as read 1, read 2, index-1, and index-2, respectively. The first 17 bases (SSR and anchor regions) in both reads were skipped using the DarkCycle option of the MiSeq system. Read 1 and 2 cannot overlap within each fragment, and the size range of the library was 400–800 bp. Thus, following to Suyama and Matsuki (2015), we treated reads 1 and 2 as independent reads. The final data output to the next step included sequences of 80 bases from both ends of 400–800 bp forward-reverse amplicons of various ISSR regions.

SNP detection

Read quality filtering was performed using Trimmomatic 0.39 (Bolger et al. 2014). Extremely short reads containing adapter sequence were filtered out. Then, to remove low quality reads, the reads were scanned with a four-base wide sliding window, and reads with an average quality per base below 15 were removed. After quality filtering, a total of 39,956,740 reads were obtained. SNP selection was performed using STACKS 2.41 (Catchen et al. 2013). First, reads were grouped to each locus using the ustacks, cstacks, sstacks, tsu2bam, and gstacks commands with the following parameters recommended by Paris et al. (2017): minimum depth option creating a stack (m) = 3, maximum distance between stacks (M) = 2, maximum mismatches between loci when building the catalog (n) = 2, and number of mismatches allowed to align secondary reads (N) = 4. After this process, the dataset of assembled loci (stacks dataset) was obtained. For most analyses (except for gene flow analysis), loci present at a rate of more than 80% among all samples were extracted using the populations command. Several options were included with this command: the minimum minor allele count over the entire dataset was set to three to exclude potentially artificial loci found in only one individual, sites showing excess heterozygosity (>0.6) were removed to filter potential heterozygotes resulting from artificial loci built from several paralogous genome regions, and the output was limited to one SNP per locus to exclude linkage disequilibrium among SNPs. After filtering, 489 SNPs were obtained.

Genetic diversity and differentiation

For each population, the expected heterozygosity (HE) and fixation index (FIS) were calculated using the populations command in STACKS. Significant deviations from Hardy–Weinberg equilibrium (HWE), as indicated by FIS deviating from zero, were tested by 1000 randomizations using FSTAT 2.9.4 (Goudet 1995). Genetic differentiation among populations was evaluated by calculating global/pairwise FST values (Weir and Cockerham 1984) using GenAlEx 6.5 (Peakall and Smouse 2012). To understand the spatial trend of the genetic structure, a Mantel correlogram displaying correlations between FST and geographic distance for each of the ten distance classes was constructed with 9999 permutations using the package VEGAN in R.

Population structure

The population structure was examined using STRUCTURE 2.3.4 (Pritchard et al. 2000), which implements a Bayesian clustering method using multi-locus allele frequency data. The STRUCTURE settings were the admixture and allele frequency correlated model with previous sampling location information (LOCPRIOR; Hubisz et al. 2009). The algorithm was run 20 times for each K from 1 to 14 with a burn-in of 20,000 followed by 30,000 Markov chain Monte Carlo (MCMC) replicates. The program CLUMPAK (Kopelman et al. 2015) was used to compile the results of the STRUCTURE analysis for each K. STRUCTURE HARVESTER (Earl and vonHoldt 2012) was employed to calculate the probability of the data for each K (LnP(D); Pritchard et al. 2000), the corresponding standard deviation, and Evanno’s delta K (Evanno et al. 2005).

Association of genetic variation with geographic distance and water temperature

To evaluate the effects of space and water temperature on genetic differentiation, the independent correlations of FST with water temperature differences and geographic distance were calculated by Mantel tests with 9999 permutations. To identify confounding effects between space and water temperature if present, multiple regression on distance matrices (MRM; Lichstein 2007) was used, with “water temperature differences + geographic distance” as the explanatory variable. Mantel tests and MRM were performed with the package ECODIST in R. These analyses were conducted at three scales: (a) using all 20 populations (waterway distance 0.3–70.6 km), (b) using upstream 12 populations (Pop1-12; waterway distance 0.3–24.1 km), and (c) using structured 9 populations (Pop4-12; populations assigned as one large cluster in STRUCTURE). We used scale (b) to account for the mobility of Cottus and to exclude the effects of a long spatial gap and a dam, and scale (c) was used to prevent biases when including sites displaying high genetic divergence (‘outliers’; Koizumi et al. 2006).

Some studies have considered it inappropriate to apply the Mantel test to the raw vector data and to take the difference to create a matrix (Legendre and Legendre 2012). Hence, we also used distance-based RDA, which is a constrained ordination approach that can use environmental variables as vector data and has higher power than the Mantel test (Harmon and Glor 2010; Legendre et al. 2015). We used genetic variation as the response variable and water temperature and spatial variables derived from geographic distance as explanatory variables. For genetic variation, a principal coordinate analysis of the pairwise FST matrix was conducted, and all axes were used as response variables. Spatial predictors were generated as a set of distance-based Moran’s eigenvector maps (MEMs; Griffith and Peres-Neto 2006), vectors that capture broad- to small-scale spatial structures. MEMs were generated from the geographic distance matrix using the package ADESPATIAL 0.3.8 (Dray et al. 2020) in R. To identify meaningful MEM predictors, forward selection (999 permutations) was used for generated MEM dataset. In each RDA with the explanatory variables of water temperature and geographic MEM variables, adjusted R-square values (Radj2), which penalize the increase in explanatory power due to an increase in the number of explanatory variables, were calculated. RDAs were performed on three scales, as in the Mantel tests and MRM analyses. RDAs including forward selections were performed with the package VEGAN in R.

Contemporary gene flow

To infer gene flow patterns, BA3SNP (Mussmann et al. 2019), a reconstructed version of BayesAss (Wilson and Rannala 2003) that estimates dual-direction pairwise migration rates (here migration is successful dispersal that contributes to gene flow) over a few generations, was used for the 20 populations. To improve the convergence of BayesAss, another SNP filtering criterion that increases the variance among populations was used. From the stacks dataset, loci present at a rate of more than 40% among at least two populations were extracted with a minimum minor allele count of three. As in the other analyses, sites showing excess heterozygosity (>0.6) were removed, and the output was limited to one SNP per locus. Then, 1291 SNPs were obtained. We set 20,000,000 MCMC iterations, including a burn-in period of 10,000,000. Following the proposal of Meirmans (2014), we performed ten independent runs with different seed values and chose the run with the lowest Bayesian deviance to overcome poor MCMC sampling. The inferred migration rates with 95% credible intervals that did not include zero were regarded as significant. We then defined the net immigration rate into a given population A, from another population B, as the [estimated value of B- >A minus that of A- >B]. Then, the mean net immigration rate for each population (the mean net immigration rate for a given population from all of its paired net immigration rate values between all other populations) was calculated to measure the degree to which a population is a donor or a recipient of migrants (Hänfling and Weetman 2006; Sexton et al. 2016). To investigate whether water temperature affects the source-sink structure independent of upstream-downstream dispersal, we constructed a linear model with a response variable of the mean net immigration rate and explanatory variables as maximum water temperature and elevation. The Pearson’s correlation between water temperature and elevation was −0.39; thus, no multicollinearity was considered. The independent effects of water temperature and elevation were assessed by the partial regression coefficient (β) and standardized partial regression coefficient (std β). Likelihood ratio tests were used to determine the significance of explanatory variables. This analysis was conducted using the package LMTEST (Zeileis and Hothorn 2002) in R.

Results

Heterogeneity of water temperature

Summer water temperatures varied within the watershed, and the maximum water temperature ranged from 12.5 to 23.5 °C. The correlation between geographic distance and water temperature difference was low (Mantel r = 0.08, p = 0.45). In the Mantel correlogram for water temperature, r values were not significant in all distance classes (Fig. 2A), indicating little spatial autocorrelation.

Fig. 2: Mantel correlograms showing spatial autocorrelation of water temperature and genetic data.
figure 2

Filled circles indicate significant correlations from a null model of spatial structure (*p < 0.05, **p < 0.01). A water temperature; B genetic variation.

Genetic diversity, divergence, and population structure

FIS did not deviate significantly from zero except in one population (Pop20), suggesting that HWE could be assumed in most populations. The HE was similar across the watershed, but one population located in the upstream section of a check dam (Pop13) displayed a slightly lower value (Table 1). Across the 20 populations, the FST value was 0.038, indicating weak population differentiation, and pairwise FST values ranged from 0.001 to 0.127 (Table S1). In a Mantel correlogram, r values were significantly positive in up to the second distance classes (~15 km) and reached zero at a distance of 27 km (Fig. 2B).

In the STRUCTURE analysis, the probability of the data (LnP(D)) increased progressively with each K, and delta K was highest at K = 4 (Fig. S1). Distinct clusters corresponding to the geographic structure were detected up to K = 6 (Fig. 1). For K = 2, the genomes of the individuals in Pop13 were grouped into a single unique cluster. For K = 4, Pop1-3 and Pop16-19 were grouped as additional clusters. Pop14,15,20 was inferred to be an admixed cluster when K = 4 and a distinct unique cluster when K = 6. Pop3 also formed a unique cluster at K = 6. Pop4-12 was assigned to almost one cluster at even large K.

Effects of space and water temperature on genetic differentiation

Mantel tests captured significant correlations between pairwise FST and waterway geographic distance (Mantel r = 0.27, p < 0.05; Fig. 3), indicating a significant effect of IBD. Conversely, no correlation between FST and water temperature differences was detected (Mantel r = −0.10, p = 0.56). The MRM with these two matrices as an explanatory variable did not show a significant relationship, and the R2 value was not very different from that of geographic distance alone. The same pattern was obtained at different spatial scales (Table 2). The forward selection of the MEM variables identified six significant predictors and the RDA with selective MEMs was significant with an Radj2 value of 0.61 (Fig. S2). Even in RDA, water temperature did not explain the genetic variation at all.

Fig. 3: Relationship of genetic differentiation with water temperature and geographic distance.
figure 3

Scatter plots of pairwise FST versus water temperature differences (A) and geographic distance (B) are displayed.

Table 2 Summary of Mantel tests, MRM, and RDAs on three spatial scales for examining the effect of space and water temperature on genetic differentiation (FST).

Gene flow

A total of 11 significant migration rates were detected (Tables 3 and S2). All of them were migration from relatively low-temperature sites (mean 14.3 °C [SD 2.90]) to high-temperature sites (mean 20.0 °C [SD 1.84]). Most of them were upstream to downstream direction, but they also included migration from downstream to upstream direction (Pop4 to 5 and Pop12 to 11). The linear model indicated that water temperature had a significant effect on the mean net immigration rates (std β = 0.542; p < 0.05) but elevation did not (std β = −0.077; p = 0.70; Table 4). Populations with lower water temperatures displayed lower mean net immigration rates, indicating that they are the source populations (Fig. 4).

Table 3 Significant migration estimated by BA3SNP (BayesAss).
Table 4 The linear model explaining the mean net immigration rate.
Fig. 4: Variation in the mean net immigration rate with maximum water temperature.
figure 4

Populations with negative net immigration (i.e., greater emigration) indicate relative source populations, whereas populations with positive values (i.e., greater immigration) indicate relative sink populations. The dashed lines indicate the 95% confidence interval.

Discussion

Genetic differentiation

In riverscapes, lower genetic diversity upstream than downstream has been recognized as a typical pattern (Blanchet et al. 2020; Thomaz et al. 2016). However, in this study, genetic diversity was similar regardless of longitudinal position (Table 1). The relative suitability of the upstream environment for the sculpin and the source-sink gene flow that is not only in the upstream-downstream direction (Table 4) may have suppressed the decrease in genetic diversity of upstream populations.

Among the entire sampling sites, a low level of genetic differentiation was observed (FST = 0.038). A study of C. nozawae in Tohoku district displayed high genetic differentiation within one watershed (Ito et al. 2018), and previous studies investigating other fluvial Cottus species within 30–50 km2 scale have often reported high FST values (Hänfling and Weetman 2006; Junker et al. 2012; Ruppert et al. 2017; but see Peacock et al. 2016). In the present study, the degree of genetic differentiation was somewhat lower than those observed in these previous studies. Other than gene flow, the large effective population size and consequent low genetic drift from the common ancestor may result in low genetic differentiation. Nevertheless, high-resolution genome-wide SNP analysis allowed genetic differentiation within the watershed to be clearly detected.

Distinct genetic structure along the geographic structure was inferred within the watershed (Fig. 1). There was a significant correlation between genetic differentiation and geographic distance, but summer water temperature had no effect on the strength of genetic differentiation (Fig. 3). Thus, the possibility that differences in water temperature lead to genetic differentiation via local adaptation is low.

Only a handful of studies have investigated the associations between the environment and the genetic structure of sculpin. Ruppert et al. (2017) analyzed the relationship between genetic and elevation differences and displayed a pattern of IBE along elevation. Another study demonstrated genetic divergence between different habitat types (Dennenmoser et al. 2014). However, most studies have been unable to eliminate the interaction of environment and geographic distance. It is not clear whether populations in these different regions and species actually have IBE patterns, but the IBE of sculpin may not be very common in wild populations.

There may be some factors related to genetic differentiation that are not clear from this study (residuals). For example, the factors promoting the genetic divergence between Pop1-3 and Pop4-12 detected in STRUCTURE analysis are still unknown. Since we assume the IBE pattern as a hypothesis, we used only an environmental variable of each site, but the factor of genetic differentiation may also lie on the pathway between populations (e.g., streambed slope). To identify those factors contributing to the residuals of genetic differentiation, analyses using environmental data across the entire watershed (e.g., White et al. 2020) and denser sampling are necessary.

Gene flow

Another important insight into population structure was obtained via gene flow analysis. Almost all significant migration was detected in population pairs within ~15 km, the range where the spatial correlation of genetic variation was significant in the Mantel correlogram (Fig. 2B). Most of the 11 cases of significant migration were downstream direction (Tables 3 and S2), which is consistent with several previous studies about fluvial sculpins that have revealed the source-sink structure from upstream to downstream sites (Hänfling and Weetman 2006; Junker et al. 2012). However, the present study also detected migrations toward tributaries at higher elevations, from lower-temperature tributaries to higher-temperature tributaries. Of course, we should be cautious about the accuracy of BayesAss estimates (Meirmans 2014). However, the estimated pattern displayed a clear trend, and it could be concluded that suitable habitats with low summer water temperatures behave as “source” habitats in the watershed. The effect of summer water temperature on the mean net immigration rate was larger than that of elevation (Table 4) and populations with lower temperatures had a greater degree of individual supply into other populations (Fig. 4).

Although BayesAss estimates migration rates over a short period, the estimates should reflect the gene flow pattern in the studied ecosystem. In this study, significant migration was estimated between populations with different summer water temperatures at close locations (Fig. 1; Table 3). This finding implies that gene flow in the direction of across dissimilar environments is high. This type of gene flow is referred to as ‘counter-gradient gene flow’ (Sexton et al. 2014). The IBE pattern is detected when gene flow among similar environments is high; therefore, the counter-gradient gene flow is the opposite pattern to IBE and suppresses IBE. The ecological processes driving this gene flow pattern could not be identified in this study, but it may at least be related to differences in population density documented by Suzuki et al. (2021) in the same watershed as the present study (i.e., high density in low-temperature streams).

Goto (1998) used mark-recapture methods and reported that most mature individuals stay at the same site for more than 1 year, indicating almost no migration. While the mark-recapture method addresses only adult movement, genetic approaches also reflect movement during early life stages that is not easily captured by mark-recapture approaches (Lamphere and Blum 2012). If passive gene flow (i.e., the transportation of juveniles or eggs) is the driving force of gene flow, some high migration rates obtained in our study do not conflict with the results of mark-recapture methods. Furthermore, some previous studies using genetic methods on Cottus reported similar migration rates to those obtained in the present study and discrepancies between genetic and mark-recapture estimates (Junker et al. 2012; Lamphere and Blum 2012), suggesting that the detected degree of gene flow is reasonable. Still, it remains to be seen whether the estimated degree of gene flow is constantly observed. For instance, flooding is a factor that increases fish dispersal (Blondel et al. 2021; Natsumeda 2003), and it should be noted that we cannot rule out the possibility that recent large disturbances (e.g., relative large typhoon in 2016) may have temporarily promoted dispersal. To make the findings on gene flow more general, more research under different conditions and other evidence such as differences in turnover rates are needed.

Implications for future management

Because sculpins are not commercially or recreationally important fish, they are impacted relatively little by human translocations (Lucek et al. 2018). Thus, sculpin genes have a footprint of a long-term past history, and the diversity of each watershed should be protected. Within the watershed, while each tributary is characterized by a unique genetic composition, loose connectivity via gene flow between tributaries is maintained. Indeed, fragmented populations such as Pop13 show lower genetic diversity. In particular, populations in streams with less suitable environments may be maintained by immigration from low-temperature streams that are more suitable habitats. Therefore, under ongoing climate change, (i) the preservation of low-temperature streams and (ii) the prevention of making unpassable barriers between source and sink populations would be the key measures for conserving sculpin populations. Considering that cold-water fish other than benthic fish have been well studied and revealed to show high mobility between tributaries, these measures will be effective for overall cold-water fish conservation. Groundwater and its associated factors such as watershed geology are useful for identifying low-water-temperature streams and conserving possible “source” populations.

Local conservation measures can also be proposed. If the distance at which the spatial coefficient becomes nonsignificant is considered to represent the appropriate management unit size (Diniz-Filho and De Campos Telles 2002), stream networks should not be subdivided into stretches of <15 km to maintain genetic resources. However, even if a large continuous habitat is maintained, then sections consisting solely of high-temperature streams will lead to local population declines due to the absence of individual supply sources. Thus, care must be taken whenever rivers are fragmented. In the Sorachi River, many C. nozawae inhabit high density (Suzuki et al. 2021), and they are not considered to be in danger of immediate extinction. However, maintaining the abundance of general species is very important to avoid ecosystem breakdown (Baker et al. 2019). In some regions in Japan (Tohoku district), C. nozawae populations are in danger of local extinction (Ministry of the Environment Government of Japan 2020). Therefore, regardless of how large the population is, maintaining continuous stream networks and the connectivity to low-temperature sections in the networks is needed for future population persistence.

Implications for landscape/riverscape genetics

Decoupling distance and environment in a locality is a strong strategy for rigorously evaluating competing hypotheses (Gilbert and Lechowicz 2004). Through this approach, we revealed a pattern that has not been proposed for the examined genus. Freshwater populations have been identified as potentially fruitful targets for the application of a landscape genetic approach to delineate population structure (Kanno et al. 2011), and we focused on groundwater to reveal the association of water temperature with the genetic structure of cold-water fish populations.

In this study, the possibility of a source-sink structure caused by counter-gradient gene flow was suggested instead of IBE. Many source-sink structures from core populations to edge (peripheral) populations have been studied. However, local gene flow patterns are also important for discussing population structure in changing environments. The IBE pattern is easily detected when spatial-environmental correlation is present, but it is important to scrutinize whether it truly exists. Although the local source-sink pattern provides notable insights for conservation and evolution (Iles et al. 2018), a previous meta-analysis showed that counter-gradient gene flows are much less frequently reported than IBE (Sexton et al. 2014). Counter-gradient gene flow is inherently exclusive to IBE and may tend to be neglected. In population genetic studies, the results often differ from the hypotheses (Myers et al. 2019). As detecting the relationship between the environment and genetic differentiation is still challenging, we hope that the role of the environment in genetic divergence will be better understood in many studies with the help of sampling strategies.