Most population genetic studies in streams infer long-term patterns of gene flow by calculating fixation indices (e.g., FST) among sampled populations. In more-recent analytical methods, the need to assign individuals to populations a priori (clustering algorithms) is relaxed, and spatial autocorrelation analysis of allele frequencies (SA) is used to infer finer-scale and potentially short-term dispersal distances. We applied multiple methods to study the population genetic structure of the riverine caddisfly Stenopsyche marmorata (Trichoptera:Stenopsychidae) from 4 adjacent catchments in northeastern Japan. We genotyped larval individuals (N = 532) from 30 sites at 8 polymorphic microsatellite loci. Fixation indices suggested low levels of genetic differentiation among populations (global FST = 0.062, p < 0.01), and significant isolation-by-distance (IBD) indicated populations were in drift-migration equilibrium. Bayesian clustering separated S. marmorata into distinct upland (>250 m asl) and lowland populations, with different FST values (upland FST = 0.048, p < 0.01; lowland FST = 0.029, p < 0.01) and significant IBD only among upland populations. Allele frequencies were significantly positively autocorrelated (Moran's I > 0, p < 0.05) at distances up to 18 km along streams and up to 12 km across terrestrial habitat. These values were similar to directly observed flight distance in a single generation for this species in the field. We conclude that the multiple-method approach revealed: 1) unexpected population subdivision between upland and lowland areas that may result from local adaptation, differences in phenology, and historical colonization by multiple lineages; and 2) fine-scale estimates of dispersal that match direct observations of flight and suggest gene flow is more pronounced along water courses in this species.
The extent of dispersal within and among suitable habitats has important consequences for the ecology, evolution, and conservation of freshwater organisms (Malmqvist 2002, Bohonak and Jenkins 2003, Hughes et al. 2009). Aquatic insects are often the predominant component of the benthic fauna in rivers and streams, and their dispersal is thought to take place largely through larval drift and adult flight. Larval drift is longitudinal and restricted to the water course, whereas adult flight may be longitudinal or lateral, with adults dispersing through terrestrial habitats to other streams. Svensson (1974) studied dispersal in one of the major groups of stream insects, the caddisflies (Trichoptera). They reported species-specific patterns that included flight in upstream, downstream, and random directions, including into surrounding terrestrial areas.Morerecent studies of caddisfly dispersal also have documented lateral movement of adults between streams of the same catchment and even across major catchment boundaries (Kovats et al. 1996, Collier and Smith 1998, Petersen et al. 2004).
Population genetic markers have been used to study caddisfly gene flow, where genetic differentiation is used as an indicator of long-term patterns of dispersal among populations. Most investigators estimated FST among already defined populations, where significant FST indicates some constraint on gene flow. Such methods are thought to integrate across multiple generations and incorporate rare or long-distance events that may not be captured by field sampling (Bilton et al. 2001, Bohonak and Jenkins 2003, Smith and Green 2005, Zimmermann et al. 2011). They also require only a single visit to the field and can be done on any life stage. Authors of several population ge- netic studies have reported evidence for gene flow within and among streams and catchments (Wilcock et al. 2003, Schultheis and Hughes 2005, Hughes 2007), and some reported relatively high levels of differentiation, results suggesting that dispersal among catchments varies according to species' flight behavior and ability (Miller et al. 2002, Watanabe et al. 2008). This approach provides an overall view of historical gene flow across long time scales, but may not give an accurate prediction of contemporary dispersal (Monaghan et al. 2001).
Alternatives to traditional population-based FST approaches include cluster analyses and analyses of statistical autocorrelation (SA) of allele frequencies. Bayesian clustering algorithms (e.g., Pritchard et al. 2000, Evanno et al. 2005) relax the need to assign individuals to populations prior to analysis, thereby providing a more objective means to identify populations before using F-statistics. SA analysis measures covariance of allele frequencies at specified distance classes (Sokal and Oden 1978a, b, Hardy and Vekemans 1999, Epperson 2005). SA has been used to estimate the spatial scale of gene flow in the short term, i.e., contemporary dispersal, by a number of taxa in various ecosystems. In simulation studies, distance classes showing positive autocorrelation were almost identical to the maximum distance of gene flow per generation (Sokal and Wartenberg 1983, Epperson 1995). Empirical studies support these findings for several freshwater organisms, i.e., the damselfly Coenagrion mercuriale (Watts et al. 2004), the salmonids Salvelinus fontinalis, Salvelinus alpinus, and Salmo salar (Gomez-Uchida et al. 2009, Kanno et al. 2011), and the salamander Plethodon cinereus (Liebgold et al. 2011). We are not aware of any studies in which the method was applied to a lotic insect, but SA analysis of distance classes in rivers and across the landscape may provide some insight into the relative propensity for dispersal along rivers and may provide dispersal distances comparable to those measured directly in cases where dispersal is difficult to measure directly.
We examined the population genetic structure and dispersal patterns of the riverine caddisfly Stenopsyche marmorata (Trichoptera:Stenopsychidae) from 4 adjacent catchments in northeastern Japan. The species is widespread and abundant in Japan, and its flight behavior (Nishimura 1966, 1967, 1981) and life cycle (Nishimura 1966, Gose 1970, Aoya and Yokoyama 1987, Kocharina 1999) have been well studied. Population genetic studies suggest the species is neither panmictic nor highly restricted, showing drift-migration equilibrium (i.e., significant isolation by distance [IBD]) at a scale of 33 km (based on analysis of random amplified polymorphic deoxyribonucleic acid [RAPD]; Watanabe et al. 2010) and reduced gene flow across 2 large reservoirs (based on RAPD; Watanabe and Omura 2007). We applied multiple methods to highly variable microsatellite loci to study population structure and compared results from traditional FST approaches with those obtained with clustering and SA analysis. We then compared these estimates of gene flow with published data on direct observations of flight in this species. We examined geographic (Euclidian) and riverine distances separately to assess whether dispersal was more pronounced along stream corridors or across the landscape.
METHODS
Study area and sampling
We conducted our study in 4 adjacent catchments with headwaters in the Ou mountains of Miyagi prefecture in northeastern Japan (Fig. 1). Streams in the area are characterized by high environmental heterogeneity along short and steep corridors. The region is largely mountainous or hilly and forested with beech (Fagus) and cedar (Cryptomeria), although catchments include agricultural areas (13% of total study area, primarily rice paddy) and residential and commercial areas (11%). We sampled 30 study sites. The Natori River (20 sites) has the largest catchment area (939 km2) and total stream length (55 km) and includes the Hirose River subcatchment (8 of the 20 Natori sites). Below the confluence, the water is brackish and S. marmorata larvae do not occur. The Nanakita River (7 sites) is the next largest catchment (229 km2). The Masuda (1 site) and Gokenbori (2 sites) rivers form relatively small catchments and flow from low hills. Four large dams, a large natural lake, a natural waterfall (Fig. 1), and ∼180 small storage reservoirs occur in the area (Miyagi Prefectural Government 2008). We conducted field surveys in July (summer) and November (autumn) 2006. At each site, we collected larvae using a kick net (250-µ;m mesh) along a 300-m stream reach. We preserved all specimens in 99% ethanol in the field, returned to them to the laboratory, and identified them to species with the aid of a stereomicroscope.
Study animal
Stenopsyche marmorata has been the subject of a number of ecological studies. Nishimura (1966, 1967, 1981) studied adult flight behavior of S. marmorata in the Yagi River (Maruyama Basin, Hyogo Prefecture, southwestern Honshu, Japan). These studies yielded quantitative estimates of flight distance that provide a benchmark for comparative estimates of dispersal. In 2 mark-recapture studies, Nishimura (1967, 1981) used fluorescent dye to mark flying adults, and set a vinyl curtain trap (3.6 m2) across the stream. He collected females 3.8 km upstream after mating but found none at a 2nd collection point 4.7 km upstream. By calculating female flight speed (2.81 m/s) based on the time it took individuals to fly between 2 points (20.5 m apart), the duration and timing of active flight (34 min in the evenings, with air temperature > 10°C, wind velocity < 5.4 m/s, and illumination < 400 lux), and the number of days individuals lived (marked individuals were captured after 2 d, but not after 3 d), Nishimura (1981) estimated that gravid females could fly up to 12 km in a generation. He also measured male flight speed (2.18 m/s) but did not calculate their flight distance. In 2 studies examining flight along the water course, Nishimura (1967, 1981) reported that 70 to 90% of adults flew upstream (30- 10% flew downstream). In a study based on collections with 1 light trap set 10 m from the stream, Nishimura (1966) reported individuals of both sexes at the trap, suggesting at least some propensity for lateral dispersal.
Figure 1.
Map of the study area indicating rivers and sampling sites. Sites are indicated by rectangles in which vertical lines indicate each individual's proportional assignment to the upland (black) or lowland (white) population (K-pop) based on the hierarchical Bayesian clustering analysis.
Larval S. marmorata are net-spinning filter-feeders that capture particulate organic matter and feed on periphyton (Nishimura 1966, Doi et al. 2007, Shin et al. 2013). Gose (1970) examined 3 populations in Japan and found that S. marmorata was univoltine in colder areas and bivoltine in warmer areas. Authors of more recent studies reported similar patterns: univoltine in colder areas (Siberia; Kocharina 1999) and bivoltine in warmer areas (Japan; Aoya and Yokoyama 1987). Aoya and Yokoyama (1987) reported that S. marmorata emerged from late May to early October with 2 distinct emergence peaks in the Su river (Mogami basin, Yamagata Prefecture, Japan), which occurs at a latitude similar to that of our study area. No reports of the emergence period have been made for our study area.
Microsatellite genotyping
We genotyped 532 individuals (mean = 17.7/site; range 11–20). We removed the digestive tract and epidermis, and isolated deoxyribonucleic acid (DNA) from the remaining tissue using DNeasy 96 Blood and Tissue Kits (Qiagen, Tokyo, Japan). We quantified DNA concentration with a spectrometer (NanoDrop-1000; Thermo Fisher Scientific, Waltham, UK). We amplified 8 microsatellite loci that we developed previously for S. marmorata (Molecular Ecology Resources Primer Development Consortium 2009) using 4 multiplex polymerase chain reactions (PCRs) (Appendix S1). We did PCR with a Veriti 96-well thermal cycler (Applied Biosystems, Foster City, California) in 10-µ;L reaction volumes with 10 ng template DNA, 0.25 U Taq DNA polymerase (TaKaRa, Tokyo, Japan), varying concentrations of primers (Appendix S1), 1 × PCR buffer (TaKaRa), 3 mM MgCl2, 0.4 mM dNTP (TaKaRa), and 1 or 2% tween20 (Appendix S1). PCR cycling conditions were 2 min at 95°C, followed by 35 cycles of 30 s at 95°C, 30 s at 47 or 50°C, 1 min at 72°C, and a final extension step of 5 min at 72°C. We analyzed fragment sizes on an ABI 3100 automated sequencer (Applied Biosystems) with Peak Scanner (version 1.0; Applied Biosystems). We estimated fragment sizes with a GeneScan 500 ROX size standard (Applied Biosystems). We confirmed the quality of all fragments with MICROCHECKER (version 2.2.3; Van Oosterhout et al. 2004).
Data analyses
We calculated the number of alleles (Na), observed heterozygosity (Ho), and expected heterozygosity (He) for each population with GENALEX (version 6.4; Peakall and Smouse 2006). We calculated allelic richness (Ra) with FSTAT (version 2.9.2; Goudet 1995). We tested for linkage disequilibrium (LD) and deviation from Hardy- Weinberg equilibrium (HWE) for each population and locus with GENEPOP (version 4.2; Raymond and Rousset 1995) with a likelihood test (settings: dememorization 1000, batches 100, iteration per batch 1000). We adjusted significance levels with Bonferroni correction. We ran a hierarchical analysis of molecular variance (AMOVA) in ARLEQUIN (version 3.1; Excoffier et al. 2005). The analysis provided estimates of genetic differentiation at 3 hierarchical levels: 1) among basins, 2) among sites within basins, and 3) within sites.
We calculated global FST for all populations in ARLEQUIN with 1000 permutations and pairwise FST with GENALEX (version 6.4; Peakall and Smouse 2006). The presence of null alleles can lead to incorrect estimates of FST by reducing the genetic diversity within populations (Slatkin 1995). Therefore, we estimated the null allele frequencies and global FST with and without the estimated null allele (ENA) correction with FreeNA (Chapuis and Estoup 2007) and the EM algorithm (Dempster et al. 1977) with 1000 permutations. For the test of IBD, we examined correlations of pairwise FST with geographical (Euclidean) distance and riverine distance (distance along a watercourse) between sites. We used Mantel tests (1- sided, 1000 randomizations) to test for significant positive relationships between pairwise FST and distance using SPAGeDi (version 1.3; Hardy and Vekemans 2002).
We delineated subpopulations with a model-based Bayesian clustering analysis of all 532 genotyped individuals as implemented in STRUCTURE (version 2.3.2; Pritchard et al. 2000). The correlated allele-frequency model (Falush et al. 2003) infers population structure where there are K subpopulations (K is unknown a priori). Given the numbers of candidate Ks, the model estimates loglikelihood (lnP[K]) for each K and an optimal K can be chosen based on the highest standardized 2nd-order rate of change (ΔK) of lnP(K) (Evanno et al. 2005). We performed 20 runs of 100,000 iterations with a burn-in of 50,000 for each candidate K ranging from 1 to 15 using the admixture model. We used a uniform prior for α, the parameter representing the degree of admixture, with a maximum of 10.0, and set Alphapropsd to 0.05. Lambda, the parameter representing the correlation in the parental allele frequencies, was estimated in a preliminary run using K = 1. The prior FST was set to the default value (mean ± SE, 0.01 ± 0.05), which was within the range of our observed global FST (see Results). Bayesian clustering identified 2 distinct populations (hereafter referred to as K-pops, see Results). Therefore, we retested for IBD and recalculated global FST and pairwise FST (as above) separately for each K-pop.
We analyzed genetic SA of alleles separately for geographic and riverine distance. We used 6 distance classes of 6 km each, based on the published estimates of flight distance for this species (above). Classes ranged from 0 to 36 km, and individuals within a site were considered 0 km distant. We calculated Moran's I for each distance class using SPAGeDi, where I ranges from -1 to 1, and positive values indicate that sites within a given distance class have similar genetic structure. We used jackknifing to estimate the 95% confidence intervals (CIs). Unequal sample size among spatial distance class can lead to bias of correlation coefficients, but the number of pairwise individuals in each distance class was similar in our analysis. We calculated SA first from all samples and then separately for upland and lowland K-pops.
RESULTS
All microsatellite loci were polymorphic (Na ranged from 10–41; Appendix S2), and sampling sites showed high levels of genetic diversity (Table 1). LD was not detected at any site or locus. Null alleles were detected in all loci (Appendix S3; frequencies in Appendix S2), but FST with ENA correction (0.058; 95% CI: 0.021–0.110) was nearly identical to FST without the correction (0.058; 95% CI: 0.023–0.106), so we did not exclude any loci from further analyses. Authors of several other microsatellite studies have reported no major difference between FST with and without ENA (Atickem et al. 2013, Ćosić et al. 2013, Heidinger et al. 2013).
Global FST among all sites was 0.062 ( p < 0.01), a result suggesting low levels of differentiation. Pairwise population differentiation varied by an order of magnitude (FST = 0.018–0.178) and significant IBD was found for all sites (Fig. 2A, B). The correlation with riverine distance (r = 0.36, p < 0.01) was higher than with geographic distance (r = 0.26, p < 0.01) (Fig. 2A, B). The Bayesian clustering analysis delineated 2 distinct populations within S. marmorata, 1 that was widespread in lowland sites and middle reaches (24 sites, all but 1 were < 250 m asl) and 1 that was found only in upland sites (5 sites, all > 250 m asl) (Fig. 1). We refer to these populations as lowland and upland K-pops. When tested separately for each K-pop, differentiation was higher in the upland (FST = 0.048, p < 0.01) than in the lowland (FST = 0.029, p < 0.01) populations. Most of the higher values for pairwise population differentiation occurred between sites from the upland and lowland K-pops. When examined separately, only upland populations (we examined only sites with >5 upland individuals: B02, C03, C05, D02, and D12) exhibited IBD (r = 0.60, p < 0.01, number of pairs = 10). Riverine distance was not analyzed for upland sites because they were not directly connected by flow. In lowland populations, IBD was not significant for geographic (r = 0.09, p > 0.05, number of pairs = 325) or riverine distances (r = 0.06, p > 0.05, number of pairs = 61). Significant differentiation existed among sites within basins (hierarchical AMOVA, df = 26, variance = 0.063) but not among basins (df = 3, variance = -0.002), results suggesting that variance between the 2 K-pops was driving this pattern.
SA was positively significant at 6-, 12-, and 18-km riverine distance classes, and at 6- and 12-km geographic distance classes using all pairs of individuals (Fig. 3A, B). Mean pairwise FST within a distance class was relatively low at short distances corresponding to high values of Moran's I, and increased with greater distance inversely with Moran's I (Fig. 3C, D). When we performed pairwise comparison using only lowland individuals, SA showed a similar trend, but with significance at shorter distance classes (6 km for geographic distance; 6 and 12 km for riverine distance).
Table 1.
Genetic diversity of Stenopsyche marmorata at each study site measured at 8 polymorphic microsatellite loci. N = number of individuals genotyped, Na = average number of alleles per locus, Ra = allelic richness, Ho = observed heterozygosity, He = expected heterozygosity.
DISCUSSION
Dispersal and spatial autocorrelation
Highly polymorphic markers are widely used to infer genetic variance within and among a priori populations of stream insects to provide a broad overview of gene flow over multiple generations (e.g., Wilcock et al. 2003, Rebora et al. 2005, Alp et al. 2012) and to study local effects of bottlenecks (e.g., Shama et al. 2011). We applied a traditional fixation-index approach and alternative methods that used the same data (microsatellite genotypes) but estimated gene flow on smaller spatial and temporal scales (SA) and relaxed the need to define populations a priori (clustering). Dispersal distances estimated using SA analysis of allele frequencies (∼12–18 km) agreed closely with estimates from direct observations of flight (∼12 km; Nishimura 1981). These results matched expectations from simulations (Sokal and Wartenberg 1983, Epperson 1995) and the few empirical studies (Watts et al. 2004, Beck et al. 2008, Liebgold et al. 2011, Kanno et al. 2011) in which short-distance gene flow was characterized by an initial positive autocorrelation that declined to 0 and then become negative at long distances. Negative values of Moran's I can arise from low migration rates (Hardy and Vekemans 1999) or the presence of very rare alleles (Epperson 2005). Studies of the damselfly (Odonata) Coenagrion mercuriale (Watts et al. 2004), brook trout (Salmonidae) Salvelinus fontinalis (Kanno et al. 2011), and the White-Winged Chough (Aves:Corcoracidae) Corcorax melanorhamphos (Beck et al. 2008) all found congruence between genetic SA and direct observation approaches. We are unaware of any previous studies of stream invertebrates that compared SA estimates with field observations. One limitation of our approach is a lack of suitable markers for stream insects. For example, only 5 other Trichoptera species have microsatellite markers available (search of nucleotide database of the National Center for Biotechnology Information [NCBI], July 2013).
Figure 2.
Isolation by distance (IBD) in Stenopsyche marmorata calculated using geographic (Euclidian) (A) and riverine (B) distance. Solid lines indicate significant correlation (Mantel test) for all data (panel A: r = 0.26, p < 0.05; panel B: r = 0.36, p < 0.05) and broken lines indicate significant correlation for upland populations (panel A: r = 0.60, p < 0.05). Upland populations are not directly connected by flow, thus squares are not present in panel B.
The SA approach has some limitations that require consideration. Treating a positive (and significant) autocorrelation coefficient at a given geographic distance as evidence for dispersal requires the assumption that the primary evolutionary process determining genetic structure is neutral rather than under selection. Here, significant IBD in S. marmorata suggested that populations are in genetic drift-migration equilibrium at the geographic scale examined (see also Watanabe et al. 2010) and, therefore, that neutral processes are most important in determining genetic structure. Another important consideration is that the number of samples for a given distance class must be equally distributed because the statistical power can be greater for distance classes with more samples (Vekemans and Hardy 2004). One way to overcome this problem is to set the distance classes a posteriori (Hardy and Vekemans 2002). We set distance categories a priori because we had a priori expectations of dispersal distance that we wanted to test based on field observations (6–12 km/generation; Nishimura1967). Nonetheless, we think that statistical power was not problematic because Moran's I was significant at all distance classes and the number of pairwise individuals did not vary widely except for the farthest distance class. Critically, I was positive, decreased to 0, and then became negative at greater distances, as has been observed by others (Gomez-Uchida et al. 2009, Kanno et al. 2011).
The significant positive SA at greater distance classes along rivers suggests S. marmorata has a greater tendency for riverine dispersal, but much evidence implies the existence of lateral dispersal as well: clear separation of upland and lowland K-pops in Bayesian clustering analysis, low differentiation among headwater populations, significant SA with geographical distance, IBD, and a relatively low global FST. Headwater sites were separated by geographic distances of 2.9–12.7 km and, therefore, fell within the range of dispersal distance estimated from SA. Authors of the few existing studies of lateral dispersal in Trichoptera (e.g., Svensson 1974, Kovats et al. 1996, Collier and Smith 1998, Elliot 2003) also found evidence for lateral migration. In summary, our findings support field observations that dispersal is primarily along water courses, but is not solely restricted to them. Miller et al. (2002) also found overall agreement between molecular inference and direct observation in 2 Trichoptera species. Species caught regularly in traps 100 m from the stream had lower overall genetic differentiation (i.e., more gene flow than species never found in the traps).
Figure 3.
Spatial autocorrelation (Moran's I ) (A, B) and mean pairwise differentiation (FST) (C, D) at 6-km distance classes based on geographic (A, C) and riverine (B, D) distance. N = number of pairwise comparisons in a distance class. Broken lines indicate the permutated 95% CI, and error bars indicate jackknifed 95% CI. All values of Moran's I are significantly different from 0 (p < 0.01).
Genetic differentiation between upland and lowland K-pops
The clear separation into upland and lowland K-pops by Bayesian clustering analysis was unexpected. Genetic differentiation across large reservoirs has been reported in this species (Watanabe and Omura 2007), but the only such reservoir in our study separated site D17 from sites D15 and D16, and these populations were all predominantly composed of lowland genotypes (Fig. 1). As a result, we conclude that the large reservoirs are not the cause of the observed pattern. We also did not observe marked changes in site characteristics (i.e., channel morphology, slope, land use, riparian cover) near the threshold elevation of 250 m asl. Our study was not designed specifically to test for such a transition, but we found no evidence that the 2 populations are locally adapted to 2 sets of environmental characteristics. Unpublished population genetic data from anonymous amplified fragment length polymorphism (AFLP) loci (KW, unpublished data) and mitochondrial cytochrome oxidase I (COI) sequences (SY, unpublished data) also indicate separation of S. marmorata into upland and lowland populations, with mitochondrial distances up to 3.0% suggesting historical separation of 2 lineages. One mechanism for the separation could be that upland and lowland populations have different phenologies. The mating period for S. marmorata is restricted to only a few days. Thus, asynchronous emergence may act as a reproductive barrier between these populations. Gose (1970) and Aoya and Yokoyama (1987) found that S. marmorata was univoltine at <90°C of monthly accumulated temperature (MAT) and bivoltine at >90°C. Available temperature data for 6 of our study sites, which were provided by Water Information System (Ministry of Land, Infrastructure and Transport, Japan, searched in 2011), suggest different voltinism patterns may be present. In upland site D15, mean MAT in the year prior to our study was 87.5°C, indicative of a univoltine life cycle. Lowland sites D16, D17, C17, and D22 ranged from 98–110°C, suggesting bivoltine life histories. Upland site D11 was below a small reservoir, reached 98°C, and was the only upland site with primarily genetically lowland individuals. Although anecdotal, this finding suggests the patterns may be a consequence of temperature-driven differences in phenology. Aoya and Yokoyama (1987) concluded that S. marmorata was bivoltine based on the presence of 2 emergence peaks and the fact that eggs and adults were found together, but their study site was at 200 m, in what may be a transition area where both upland and lowland K-pops occur (similar to our site C03; Fig. 1).
Conclusions
We combined traditional FST and IBD approaches to study population genetic structure with Bayesian clustering and SA of allele frequencies. Traditional methods revealed a low degree of genetic structure throughout the study area and a significant drift-migration equilibrium, whereas SA analysis revealed more fine-scale estimates of gene flow that largely agree with direct observations of flight distance in a single generation for this species (Nishimura 1981). The clear separation of S. marmorata into upland and lowland populations was not detected with traditional analyses, but plays an important role in determining genetic structure. The cause of this separation is not known, but a working hypothesis is that it is the result of historical divergence into separate lineages with different phenologies.
ACKNOWLEDGEMENTS
This research was partially supported by the Japan Society for the Promotion of Science (JSPS) (grant numbers 24254003, 25303020, 25289172). SY was supported by a JSPS Research Fellowship for young scientists, by a JSPS postdoctoral fellowship for research abroad, and by a Leibniz Association travel grant. KWwas supported by aMarie-Curie International Incoming Fellowship (EU-FP7, Co 237026). Part of this research was carried out while MTM was at Tohoku University supported by a JSPS Invitation Fellowship (S-11074). K. Muraoka and Y. Shinotsuka graciously provided the use of laboratory facilities and provided assistance genotyping. Y. Kikuchi and Y. Hamamoto helped with field work. We sincerely appreciate the comments and suggestions of Debra Finn and 2 anonymous referees, which greatly improved the manuscript.
