Evidence for contemporary and historical gene flow between guppy populations in different watersheds, with a test for associations with adaptive traits

Abstract In dendritic river systems, gene flow is expected to occur primarily within watersheds. Yet, rare cross‐watershed transfers can also occur, whether mediated by (often historical) geological events or (often contemporary) human activities. We explored these events and their potential evolutionary consequences by analyzing patterns of neutral genetic variation (microsatellites) and adaptive phenotypic variation (male color) in wild guppies (Poecilia reticulata) distributed across two watersheds in northern Trinidad. We found the expected signatures of within‐watershed gene flow; yet we also inferred at least two instances of cross‐watershed gene flow—one in the upstream reaches and one further downstream. The upstream cross‐watershed event appears to be very recent (41 ± 13 years), suggesting dispersal via recent flooding or undocumented human‐mediated transport. The downstream cross‐watershed event appears to be considerably older (577 ± 265 years), suggesting a role for rare geological or climatological events. Alongside these strong signatures of both contemporary and historical gene flow, we found little evidence of impacts on presumably adaptive phenotypic differentiation, except perhaps in the one instance of very recent cross‐watershed gene flow. Selection in this system seems to overpower gene flow—at least on the spatiotemporal scales investigated here.

. The importance of historical gene flow for ongoing adaptation is less well understood. On the one hand, we might expect the power of selection to quickly overcome any historical legacy, such that adaptive trait divergence bears little relationship to neutral genetic marker divergence (Merilä & Crnokrak, 2001). On the other hand, some studies have suggested that populations coming from different colonization events can maintain important differences in adaptive traits despite long-term occupancy of similar environments, an effect referred to as historical contingency (Losos, Jackman, Larson, Queiroz, & Rodriguez-Schettino, 1998;Taylor & Donald McPhail, 2000;Travisano, Mongold, Bennett, & Lenski, 1995). Although numerous studies have investigated whether phenotypic variation correlates with contemporary selection or with historical contingency (Alexander, Taylor, Wu, & Breden, 2006;Hoekstra, Krenz, & Nachman, 2005;Thorpe, Malhotra, Black, Daltry, & Wuster, 1995), the focus has not been on effects of contemporary versus historical gene flow on adaptive traits. To gain insight into such effects, we combine genetic inferences about historical and contemporary gene flow in guppies from two adjacent watersheds in Trinidad, with information on an important class of adaptive traits-male color.
In keeping with the expectation that the greatest barrier to dispersal in such systems is dry land, greater genetic differences are usually found among rather than within watersheds (Barson et al., 2009;Carvalho et al., 1991;Suk & Neff, 2009;Willing et al., 2010).
However, exceptions are known wherein guppies occupying some tributaries in one watershed can show surprising genetic similarity to particular populations in other watersheds (Willing et al., 2010).
These cross-watershed affinities could reflect historical or contemporary gene flow owing to natural events, such as earthquakes or severe flooding, or human-mediated transport. The best predictors of such cross-watershed gene flow events are expected to be similar elevations and geographic proximity, except in the case of some longer-distance human-mediated transfers.

| Our focal study system
Our work focused on guppies located on the north slope of the Northern Mountain Range in Trinidad, where multiple streams run roughly parallel to each other from the mountains over a series of waterfalls to the ocean. We studied two neighboring watersheds, the Marianne and the Paria (Figure 1). Gene flow within these watersheds is expected to be relatively high, at least in the downstream direction, as their total lengths are only 10.69 km for the Marianne and 9.22 km for the Paria. However, gene flow can be reduced in the upstream direction owing to the direction of water flow and to physical features such as waterfalls (Crispo et al., 2006). Waterfalls are present throughout the entire course of the Marianne, but there is no major waterfall along the Paria watershed that could prevent upstream migration. Gene flow would seem less likely between the two watersheds, and yet still might be possible owing to their close proximity at two elevation ranges: 50-100 m (downstream area) and 550-600 m (upstream area) (Figure 1). To examine patterns of gene flow with these expectations and possibilities in mind, we analyzed guppies from several sites for variation at 10 and 42 microsatellite loci.
F I G U R E 1 Topographic map of the Marianne and Paria watersheds. Bold and circled sections of the rivers indicate potential gene flow zones between watersheds, located at different elevations To see how gene flow might relate to adaptation, we focused on male guppy color, which has a known genetic basis (Gordon, López-Sepulcre, Rumbo, & Reznick, 2017;Lindholm & Breden, 2002) and evolves in response to sexual selection favoring conspicuousness and natural selection favoring crypsis (Endler, 1980;Reznick, Bryga, & Endler, 1990;Reznick & Endler, 1982). Most obviously, populations in different predation environments show dramatic color differences that reflect adaptation to the local balance between natural and sexual selection (Endler, 1980;Kemp, Reznick, & Grether, 2008;Millar, Reznick, Kinnison, & Hendry, 2006). These color patterns often (but not always) evolve quickly (Gordon et al., 2017;Karim, Gordon, Schwartz, & Hendry, 2007;Kemp et al., 2008), and differences among populations are stable over time (Gotanda & Hendry, 2014). Large variation also exists among populations within a given predation regime, reflecting the specific types and densities of local predators (Endler, 1978;Kemp et al., 2008;Millar & Hendry, 2012;Millar et al., 2006;Weese, Gordon, Hendry, & Kinnison, 2010), canopy cover (Grether, Millie, Bryant, Reznick, & Mayea, 2001;Schwartz & Hendry, 2010), and sexual selection (Houde & Endler, 1990;. Divergence in male color also has been argued to be influenced-both positively (increased divergence) and negatively (decreased divergence)-by gene flow (Endler, 1978;Fitzpatrick et al., 2015Fitzpatrick et al., , 2017. We therefore compare patterns of gene flow with patterns of neutral genetic differentiation to infer the potential role of historical and contemporary gene flow in shaping adaptive trait divergence. In summary, our goals were to (a) investigate population genetic structure of guppies in the two watersheds, (b) infer the existence and timing of gene flow events between sites within and between watersheds, and (c) test for associations between gene flow and differences in male color. Our interpretations proceed as follows: 1. If current watershed structure is the primary determinant of gene flow, samples should cluster by watershed; and gene flow estimates should be higher within than between watersheds. Deviations from this expectation (e.g., some clustering and inferred gene flow between watersheds) would indicate cross-watershed genetic exchange.

2.
If inferred gene flow between watersheds was due to historicaland presumably rare-events, such as earthquakes or floods, estimated divergence times between sites should be older than a few centuries. Deviations from this expectation (e.g., more recent divergence) would suggest the importance of contemporary factors, such as recent natural or human-mediated gene flow.

3.
If gene flow among sites is influencing adaptation, we expect patterns of male color divergence among sites to be associated with patterns of neutral genetic divergence. Deviations from this expectation (e.g., limited or no correspondence between neutral and adaptive divergence) would inform the extent to which local selection overcomes historical and contemporary gene flow, or would indicate genetic drift.

| Fish sampling
We sampled fish along the Marianne and Paria watersheds in northern Trinidad over several years (2002-2014; average of 38 individuals per year and site, min = 18, max = 50; details in Supporting Information). At each site, we used butterfly nets to capture guppies that were then transported to our laboratory in Trinidad. The fish were euthanized with a solution of tricaine methanesulfonate (MS222) and preserved in 95% ethanol for genotypic analysis. A subset of the fish was also photographed following a standard method (details below).

| Genotypic data
Two datasets were generated using different methods implemented at different times. One dataset has fewer loci (10) but more sites (20), whereas the other dataset has more loci (42) but fewer sites (12). The term "site" refers to a specific discrete sampling location, and site numbers correspond to those reported in previous work on these watersheds (Crispo et al., 2006;Gotanda et al., 2013;Gotanda & Hendry, 2014;Millar et al., 2006;Schwartz & Hendry, 2010). We here analyze both datasets because they represent two independent efforts, with different strengths and weaknesses, to quantify genetic population structure for the same watersheds. We analyze the two datasets separately because few of the loci and only some of the sites were in common.
For the second dataset ("42loci-12sites"), DNA was extracted using the same method for 42 di-and trinucleotide microsatellite loci selected from the guppy genome (NCBI BioProject PRJNA238429).
The 42 loci were multiplexed in a single PCR, and indexing sequences were subsequently added to the PCR products using a second PCR.  Raymond & Rousset, 1995) was used to test for deviations from Hardy-Weinberg equilibrium (HWE) for each locus in each "sample" (i.e., a microsatellite dataset at a particular site in a particular year) and to test for linkage disequilibrium between loci within each sample. Lositan was used to check whether l oci were potential l y under sel ection based on an F ST outlier test (Antao, Lopes, Lopes, Beja-Pereira, & Luikart, 2008).
We used the R software with RStudio (R Core Team, 2018;RStudio Team, 2016) to calculate basic population genetic measures.
The package pegas (Paradis, 2010) and package hierfstat (Goudet, 2005) were used to calculate the number of alleles per site, as well as observed heterozygosity and gene diversities (H o and H s , respectively). The package hierfstat (Goudet, 2005) was also used to calculate pairwise F ST between samples.
structure (version 2.3.4;Pritchard, Stephens, & Donnel l y, 2000) was used to infer genetic population structure and to find the appropriate number of clusters (K) that best explain the genotypic distribution. Three iterations were run for each K, from 1 to 28 or from 1 to 19 (total number of samples from the two datasets). Burnin length and run length of the program were each set at 100,000 using the admixture model and the correlated alleles model. We used the Evanno, Regnaut, and Goudet (2005) method implemented in structure harvester to find the best K. We generated structure pl ots using the R package pophelper (Francis, 2017). These analyses included multiple years of sampling for a given site so as to help assess the temporal consistency of among-site patterns.
We estimated long-term gene flow by calculating migration rates (M = m/µ) between sites using Migrate (version 3.6; Beerl i, 2009). In cases where a dataset included multiple years from a single site, we kept-for ease of estimation-only one year per site by choosing samples with the minimum length of time between them (i.e., temporal outliers were more likely to be excluded). We used an MCMC with Bayesian inference coalescent approach that employed a Brownian model approximating a single-step mutation model and default values from the software. A mutation rate of 5 × 10 −4 was chosen because it is the mutation rate commonly used in other microsatellite fish studies (Barson et al., 2009;Lippé, Dumont, & Bernatchez, 2006). For each dataset, a first run determined F ST parameters that were then used as start values for three more runs.
The number of runs was dictated by when the mean across runs was stable.
We explored recent (over the last few generations) migration rates using Bayesass (version 3; Wilson & Rannala, 2003). For each dataset, we only kept one year per site, and we first adjusted the mixing parameters to meet acceptance rates. The burn-in period of the model was then set at 1 × 10 6 , while MCMC iterations were set at 1 × 10 7 . We ran several instances of the model with different starting seeds: Results were similar among runs and so we here report only values from the first run. Model convergence was also tested using tracer (version 1.6; Rambaut & Drummond, 2013). Values calculated with this method represent the fraction of individuals in a population that are migrants derived from another population.
We used DIYABC (version 2.1.0; Cornuet et al., 2014) to estimate divergence time between pairs of sites across watershedsthe level at which such inferences were desired. This analysis was done using the 42loci-12sites dataset, with only one year per site.
For each pair of sites, we tested a simple model of two populations having diverged t generations in the past from a common ancestral population (Supporting Information), a reasonable approximation of a discrete cross-watershed gene flow event. Our models thus simplify a complex scenario of watershed colonization with multiple sites but allow the comparison of pairs of sites across watersheds. The mutation model was left as the default in the program (mean mutation rate: 5 × 10 −4 ). We generated 1 × 10 6 simulated datasets to estimate the divergence time between each pair of sites. As guppies can have 2-3 generations per year (Magurran, 2005), we assumed a value of 2.5 generations per year.

| Phenotypic data
Differences in color among guppy populations in the studied watersheds are remarkably consistent through time (Gotanda & Hendry, 2014), and so we were able to use phenotypic data (male color) from years other than the genetic data. Specifically, the data re-analyzed here were previously published in Millar et al. (2006), wherein details are provided. In brief, we extracted color information from standardized digital pictures of male guppies. Scion Image (Scion, 2001) was used to measure body size (area, length, and depth) and each color spot (area) on the left side of the body. Each color spot was classified into a color category: orange (includes red), black, fuzzy black, yellow, blue (includes purple), green, violet-blue, bronze-green, and silver. For simplicity, these categories were then further grouped into three more inclusive categories: melanic (black and fuzzy black), carotenoid (orange and yellow), and structural (blue, green, violet-blue, bronze-green, and silver). These categories and labels are only general as, for example, the "carotenoid" colors include additional compounds influencing color (Grether, Hudon, & Millie, 1999). For the present analysis, we used-for each individual fish-the total number of spots and the relative total spot area (total spot area divided by the total body area of the fish) for each color category.
We used a MANOVA to detect differences in male color between predation regimes in the Marianne watershed. We calculated pairwise P ST as a measure of the phenotypic (color) distance between guppies at each pair of sites. Following Phillimore et al. (2008), we used the formula P ST = 2 GB ∕( 2 GB + 2 2 GW h 2 ), where 2 GB and 2 GW are the between-and within-group variance and h 2 is the heritability. Given the established very strong genetic basis for the sorts of traits measured here (Gordon et al., 2017;Karino & Haijima, 2001;Lindholm & Breden, 2002;Tripathi et al., 2009), we made the assumption that h 2 = 1, meaning that all variance is genetic. Choice of a different value for heritability would not have influenced conclusions, which are based on relative differences between various types of population pairs. Following Phillimore et al. (2008), we conducted pairwise MANOVA for all sites across watersheds using the R package stats.
Variance-covariance matrices were then summed to estimate 2 GB and 2 GW .

| Comparison of genotypic and phenotypic data
To enable direct comparisons of population structure between the genetic (from both datasets) and phenotypic (male color) data, we analyzed both types of data using discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010) implemented in the R package adegenet. This method infers individual exchangeability between sites and allows evolutionary inferences from the classification of each individual into different categories of sites (Hendry, Kaeuffer, Crispo, Peichel, & Bolnick, 2013). For each data type (genetic or phenotypic), we recorded the probability that each individual is assigned to (a) its site of origin, (b) a site from the same watershed at the same elevation (upstream vs. downstream), (c) a site from the same watershed but with a different elevation, (d) a site from the other watershed with the same elevation, and (e) a site from the other watershed with a different elevation. We then recorded "classification" as the highest assignment of an individual to its own site or in any other site, and "cross-classification" as the highest assignment to any other site apart from the site of origin .
We calculated F ST (each dataset separately) and P ST means and confidence intervals to allow comparisons and used a Mantel test in the R package vegan (Oksanen et al., 2018) to statistically compare these measures. Here, we used only F ST measures from the 10loci-20sites dataset, because insufficient overlap occurred between sites in the 42loci-12sites dataset and the male color dataset.

| RE SULTS
We start with a brief summary of the main findings and the analyses supporting them before moving to specific presentation of the specific analyses. Overall, we found strong evidence of gene flow not only within watersheds but also between them-as supported by five analyses. First, structure most strongly supported four clusters for the 10loci-20sites dataset and three clusters for the 42loci-12sites dataset, with one of the clusters in each structure analysis including sites from both watersheds. Second, DAPC for the genetic data showed that, although individuals were mostly classified to their site of origin, a reasonable number were classified into sites in the same watershed (especially at the same elevation), and some were even classified into sites in the other watershed (almost exclusively at the same elevation). Third, Migrate suggested considerable historical gene flow both within and between watersheds, with estimates between watersheds often higher than some of those within watersheds. Fourth, Bayesass suggested very recent migration events among some sites within watersheds, and even from the Paria to Marianne (especially at higher elevations). Fifth, DiyaBc inferred historical gene flow between the two watersheds at low elevations and recent gene flow between the watersheds at higher elevations.
Finally, we found that male color showed little or no correspondence with neutral genetic variation, suggesting that selection tends to erase the effects of gene flow in these particular comparisons.
Although these deviations were haphazardly distributed across loci and samples, we nevertheless searched for correlations between F IS and F ST at the level of individual loci (Waples & Allendorf, 2015). We did not find any positive relationships between these two measures (10loci-20sites: r 2 = 0.08; 42loci-12sites: r 2 = 0.01), and thus ruled out a Wahlund effect. F IS values for Pre26 were very positive compared with other loci (Pre26: F IS = 0.23; median for the other loci: F IS = 0.04), reflecting heterozygote deficiency. For this locus, Microchecker indicated potential null alleles in 11 of the 28 samples. This locus was thus excluded from further analyses. In the 42loci-12sites dataset, 55 out of 798 tests showed departures from HWE after sequential Bonferroni correction. Because these departures only constituted 6% of the tests, and were haphazardly distributed across loci and sites, we did not exclude any loci from this dataset.
In the 10loci-20sites dataset, 3 out of 1,260 tests showed evidence of linkage disequilibrium after sequential Bonferroni correction. All significant tests were for site P13, which could indicate a small effective population size (N e ), or could show admixture of different lineages for guppies at that site. In the 42loci-12sites dataset, 15 out of 8,436 tests showed evidence of linkage disequilibrium after sequential Bonferroni correction. However, physical linkage is unlikely given that the loci are known (i.e., specifically developed) to be widely distributed in the guppy genome.
The F ST outlier method implemented in Lositan detected six loci potentially under selection in the 10loci-20sites dataset, and 16 loci in the 42loci-12sites dataset. For the 10loci-20sites dataset, we did not eliminate any loci, because of low information with only four remaining loci for the analysis. For the 42loci-12sites dataset, we ran structure with and without the potentially selected loci and obtained the same results. Hence, we kept all loci for subsequent analyses.
The total number of alleles per site ranged from 34 to 141 for the 10loci-20sites dataset, and from 54 to 262 for the 42loci-12sites dataset (Supporting Information). Mean number of alleles per locus was 27.11 for the 10loci-12sites dataset and was 13.02 for the 42loci-12sites dataset. Observed heterozygosity ranged from 0.292 to 0.752 for the former and from 0.073 to 0.573 for the latter (Supporting Information). Average observed heterozygosity was higher in downstream sites (H o = 0.73 ± 0.06; 10loci-20sites dataset) than in upstream sites (H o = 0.54 ± 0.14; 10loci-20sites dataset). F ST values were higher between sites that were geographically more distant (Supporting Information).

| STRUCTURE
The most likely number of clusters was four for the 10loci-20sites dataset and three for the 42loci-12sites dataset (Figure 2  P13, P14, P16, and P18). For sites that were sampled multiple times, we found consistent patterns between years. Considerable admixture between the clusters was inferred for P1, P18, M7, and M15, and admixture increased between years for M7 in the 10loci-20sites dataset. Summarizing these patterns, sites did not cluster together exclusively by watershed but rather also according to their geographic position (upstream vs. downstream; and eastern vs. western in the Marianne). We also found moderate support for a structure of 15 clusters for the 10loci-20sites dataset and 10 clusters for the 42loci-20sites dataset (Figure 2). These clusters are much more conservative; that is, for both datasets each site often constitutes its own exclusive cluster, with the notable exception of sites located in a portion of the river that is called the "Petite Marianne" (M9, M10, and M11), which cluster together in both datasets.

| DAPC
For both genetic datasets, classification was highest to the site of origin (Figures 3 and 4), indicating that each site constitutes its own guppy population at least partially isolated from other guppy populations. Some individuals were also assigned to sites from the same watershed at the same elevation, presumably reflecting the easiest contemporary routes for ongoing gene flow. For cross-classification (i.e., assigning all individuals away from their population of origin), the highest classification was generally to sites in the same watershed at the same elevation, then to the same watershed at a different elevation or instead to the other watershed at the same elevation (Figure 4).

| BAYESASS
Both datasets suggested reasonable levels of contemporary gene flow between pairs of sites in the same watershed (Supporting Information). Some sites obviously received considerable migrants from neighboring sites; for example (10loci-20sites) from P1 to P3, P7 to P17, P13 to P12 and P14, P15 to P17, M10 to M9 and M11; and (42loci-12sites) from P1 to P18, P7 to P15 and P7, P15 to P7, M7 to M8, M8 to M7 and M10, and M10 to M7. Evidence of cross-watershed contemporary gene flow was also apparent in the 42loci-12sites dataset, with some sites apparently receiving relatively recent migrants from sites in the adjacent watershed; for example, from P7 to M3 and M4 (upstream Paria to upstream Marianne), from P18 to M7 and M8 (downstream Paria to downstream Marianne), and from M9 to P18 (downstream Marianne to downstream Paria; Supporting Information). We are not certain whether these reflect actual contemporary gene flow events or rather the continued signature of past gene flow events.

| DIYABC
Divergence time estimates between watersheds differed greatly among the various pairs of sites ( Figure 5). The shortest divergence F I G U R E 3 Scatter plot based on discriminant analysis of principal components for (a) 10loci-20sites neutral markers, (b) 42loci-12sites neutral markers, and (c) male color traits. Colors correspond to a posteriori groups defined by the DAPC analysis. Individuals are represented as dots and groups as inertia ellipses time (41 ± 13 years) was estimated between nearby sites located in the upstream reaches of the two rivers. The second shortest divergence times were estimated between the downstream Paria and the upstream Marianne (533 ± 167 years) and between the adjacent downstream reaches of the two watersheds (577 ± 265 years). The longest divergence times (2,803 ± 470 years) were estimated between M16 in the western Marianne and various sites in the Paria.

| Genetic versus phenotypic patterns
General patterns here were several. First, male color patterns significantly differed between predation regimes (MANOVA; Wilks' λ = 0.766, df = 290, p < 0.001). Second, classification in DAPC was always highest to the site of origin in all datasets (Figure 4), indicating that each site is a unique "population" to at least some extent. Third, populations differed less phenotypically than genetically at all levels, especially across watersheds (Figures 3 and 6, and Supporting Information). This outcome was mostly driven by variation in neutral genetic differentiation ( Figure 6).
Together, these results suggest that phenotypic differentiation, while present among all sites, is ultimately more "constrained" in the magnitude of divergence. Third, no correspondence was seen between patterns of neutral genetic differentiation and patterns of phenotypic differentiation

| D ISCUSS I ON
Many studies have emphasized the particular nature of connectivity in streams through ideas such as the "river continuum concept" F I G U R E 4 Ratio of the mean number of individuals classified into each category as indicated on the x-axis to the mean number expected to be classified into these categories by chance. Upper panel shows the classification for the three datasets (10loci-20sites in dark gray, 42loci-12sites in medium gray, and male color in light gray), and lower panel shows cross-classification  (Vannote, Minshall, Cummins, Sedell, & Cushing, 1980). Exchanges pies has also provided evidence for this type of genetic structure.
For example, different watersheds tend to have genetically divergent guppy populations (Alexander et al., 2006;Barson et al., 2009;Carvalho et al., 1991;Suk & Neff, 2009;Willing et al., 2010) and upstream guppy populations show reduced genetic variation consistent with rare colonization events, limitations to upstream gene flow, and possible frequent bottlenecks due to floods (Barson et al., 2009;Crispo et al., 2006;van Oosterhout et al., 2006). At a first cut, our analyses suggest much the same, with the largest among-site genetic differences occurring between watersheds, with upstream Some of these linkages can be attributed to known human-mediated introductions , whereas others are more mysterious (Suk & Neff, 2009;Willing et al., 2010). In our case, cross-watershed gene flow was found in two areas, both closely adjacent tributaries at the same elevation. This finding is reminiscent of a recent study of Atlantic salmon (Salmo salar), where fish located at the same elevation but in different rivers were more related to each other than to fish in the same river but at a different elevation (Cauwelier, Stewart, Millar, Gilbey, & Middlemas, 2018). Uniquely in our study, the two cross-watershed genetic linkages seemed to have occurred on different timescales (historical and contemporary) in different parts of the watershedsneither of which are associated with any known human-mediated introductions.

| Contemporary and historical cross-watershed gene flow
We We also found signatures of historical gene flow between the two watersheds. Such signatures have been documented for some other systems and have been attributed to rare and severe events such as ice dams, earthquakes, or volcanic activity (Burridge, Craw, & Waters, 2006Gelmond, Hippel, & Christy, 2009;Lescak et al., 2015). In Trinidad's Northern Range, we inferred historical Third, guppies are currently found in the headwaters of the Jordan River, less than 50 m of horizontal distance, with an only minor F I G U R E 6 Comparison of the F ST values for both genetic datasets and the P ST values for the male color traits, within sites located in the Marianne (all sites), within sites located in the Paria (all sites), between sites located in the upstream reached of the watersheds ("recent" gene flow event, between P7 and M3-M4), sites located in the downstream reaches of the watersheds ("old" gene flow event, between P18 and M9-M10), and finally all "other sites across" watersheds other. However, this explanation seems unlikely given the remote location of these small tributaries, and the fact that indigenous people relied mainly on fish from the ocean rather than fresh water (Newson, 1976). Second, Trinidad is located on the Caribbean tectonic plate, and major earthquakes have been reported since written history of the island (Shepard & Aspinall, 1983). Such earthquakes, violent hurricanes, or massive flooding could have led to river capture (Bishop, 1995), that is, "the transfer of part or all of a (generally well established) river's flow to another river," causing the movement of Paria guppies from the Jordan River into the Petite Marianne.
Once cross-watershed gene flow occurs, a logical question is whether that influence spreads far beyond the site of origin. Several studies have shown that experimental introductions of guppies have genetic consequences that spread downstream, including over waterfalls and into different predation environments (Becher & Magurran, 2000;Fitzpatrick et al., 2015;Fraser, Künstner, Reznick, Dreyer, & Weigel, 2015). For our non-experimental, whether natural or anthropogenic, cross-watershed transfers, we also see signa-

| Consequences for adaptive traits
We uncovered signatures of gene flow within and between riverine networks reflecting a complex combination of water flow (biased downstream), barriers (waterfalls), physical proximity, potential recent human introductions, and past geological or climatological events. To what extent has this gene flow influenced adaptive trait variation? A classic theoretical expectation would be reduced divergence in the case of very high, and especially recent, gene flow (Hendry, Day, & Taylor, 2001;Lenormand, 2002). On the other hand, some theoretical treatments suggest a potential positive role for gene flow in facilitating local adaptation (review: Garant et al., 2007).
Overall, male color was quite location-specific (Figure 4), suggesting adaptation to local conditions. Some of this variation was associated with differences in predation regime (high vs. low) within the Marianne, as described in previous analyses of this system (Gotanda & Hendry, 2014;Millar et al., 2006). Yet considerable variation was also seen between our study sites within a given predation regime (Figures 3 and 4), which previous studies have attributed to these site-specific factors such as specific predator identities and densities (Millar et al., 2006), canopy cover (Grether et al., 2001;Schwartz & Hendry, 2010), and sexual selection .
However, differences among sites in color were generally less than differences among sites in neutral markers ( Figure 6 and Supporting Information). This result likely reflects some level of convergent stabilizing selection on male color owing to constraints on the range of possible color space and the need to be attractive to females but also cryptic to predators. By contrast, neutral markers are free to diverge to an extent (mostly) unconstrained by selection, instead being limited only by time and gene flow.
Importantly, we see little evidence that the constraint imposed on divergence for male color reflects gene flow-given the overall lack of correspondence between genetic and phenotypic diver- entiation. An alternative possibility is that the environments experienced by these two populations were exceptionally similar, and thus favored similar phenotypes, although habitat data do not suggest such exceptional similarity (e.g., Millar et al., 2006). Furthermore, we cannot rule out that gene flow constrains or facilitates adaptation for other traits or in other contexts. For instance, Fitzpatrick et al. (2017) found evidence for trait-specific constraining and diverging effects in an experimental manipulation of gene flow.

| Divergence time
Divergence time between the Paria and Marianne guppies was previously estimated to be approximately 100,000 years based on mitochondrial DNA data (Fajen & Breden, 1992). Our multilocus estimates suggest much more recent connections between the two watersheds, ranging from a maximum of a few thousand years between isolated portions of the watersheds, up to contemporary gene flow between proximate portions of the watersheds at the same elevation. These results have to be tempered because we only tested very simple models of divergence and because homoplasy occurs with microsatellite markers, which could create noise in our results (Estoup, Jarne, & Cornuet, 2002). However, in light of our findings, we would still like to discuss that divergence time between other watersheds extensively studied in Trinidad might also be more recent than previously estimated, a possibility that has important implications for our understanding of adaptive evolution and early speciation in this system. For instance, the general lack of speciation in guppies is often considered surprising (Magurran, 1998) given their ancient divergence-but perhaps gene flow has been much more recent. Also, although we know through experiments that contemporary evolution is common in guppies (Reznick et al., 1990), perhaps even naturally established populations have evolved on much shorter than expected time frames.

| CON CLUS ION
Our findings are broadly consistent with previous population genetic work for riverine organisms in general, and for guppies in particular. Specifically, we confirmed within-watershed gene flow in which upstream populations are less genetically diverse and more isolated than are downstream populations. However, we also discovered levels of cross-watershed gene flow, to the extent that some populations are more closely related to populations in the adjacent watershed than they are to some populations within their own watershed.
Although surprisingly genetic similarities between the Paria and the Marianne watersheds have been previously suggested (Willing et al., 2010), our much more detailed sampling was able to infer where, when, and in what directions these genetic exchanges took place.
In one case, cross-watershed linkages were recent and, in the other case, they occurred centuries ago, suggesting different contributions from geological, climatological, or anthropogenic drivers. However, none of these gene flow patterns seemed to have any major consequence for adaptive trait variation-although our findings do not rule out effects for other traits or on smaller spatial scales. Dispersal, and thus subsequent gene flow, clearly paves the way for colonization of new environments, but it did not seem to here substantially constrain adaptation by guppies to those environments.

ACK N OWLED G M ENTS
We thank three anonymous reviewers and the editor. We thank M.