Do genomics and sex predict migration in a partially migratory salmonid fish, Oncorhynchus mykiss?

We studied partially migratory populations of O. mykiss in two tributaries of the South Fork Eel River, Elder and Fox creeks, both located within the University of California Angelo Coast Range Reserve (Mendocino County, California; Fig. 1). Elder Creek drains 16.8 km² and Fox Creek drains 2.7 km², and both are steep, shaded streams. There is a waterfall located 2 km from the mouth of Elder Creek that is passable to adult steelhead under some stream flow conditions, and steelhead have been observed spawning above this barrier (Trush 1989). Oncorhynchus mykiss represent >99% of the fish biomass in these streams. The only other fish species encountered is the occasional Pacific lamprey (Entosphenus tridentatus).

We sampled fish in Fox and Elder creeks from late July to early August in 2014–2017. Fish were collected using three-pass backpack electrofishing. We sampled approximately 20% of the pools in each stream, with sample pools distributed longitudinally from the mouth to the upper extent of fish in each stream. Sample pools were selected using spatially stratified random sampling in 2014, and the same sites were revisited in 2015–2017. At capture, fish were measured for fork length (FL, mm) and mass (0.01 g), and a small tissue sample was removed from the caudal fin for genetic analyses. Fish that exceeded 2 g and 60 mm FL were tagged with a 2 mm passive integrated transponder (PIT) tag, which allowed us to track life history expression. All fish were scanned for PIT tags prior to tagging, and any recaptures were noted and remeasured for length and mass.

We extracted DNA from caudal fin samples as described by Ali et al. (2016). We included all of the samples collected in 2014 and a subsample from 2015–2017, where every other sample pool was included in the genetic analyses. We conducted restriction site-associated DNA (RAD) capture (Rapture) following the methods and bait set described in Ali et al. (2016). Briefly, the bait set targets 500 SbfI RAD tags that are distributed across all chromosomes. While the capture enriches coverage at the targeted loci, many other RAD tags are also sequenced at lower coverage. This approach facilitates analyses that require high coverage using the targeted loci while also providing a larger amount of data for analyses that do not require high coverage (e.g., discriminant analysis of principal components (DAPC) analysis from single-read sampling; see below). Libraries were sequenced using paired-end 100-base pair (bp) (samples collected in 2014) or 150-bp (all other sample years) reads on Illumina HiSeq 2500 or HiSeq 4000 machines. Demultiplexing of sequence data was performed by requiring a perfect barcode (unique to each sample) and partial restriction site match (Ali et al. 2016). Sequences were aligned to a recent rainbow trout genome assembly (https://www.ncbi.nlm.nih.gov/assembly/GCF_002163495.1/) using the MEM algorithm of Burrows–Wheeler Aligner (Li and Durbin 2009) with default parameters. SAMtools (Li et al. 2009) was used to filter alignments for proper pairs, sort alignments, remove polymerase chain reaction (PCR) duplicates, and index binary alignment map files.

All Rapture sequencing data analyses were performed using analysis of next generation sequencing data (Korneliussen et al. 2014) with a minimum mapping quality score of 20 and a minimum base quality score of 20. To select sites (SNPs) appropriate for downstream analyses, the following steps were applied. Using all samples, major and minor alleles were inferred for sites with a high probability of being variable (SNP p value < 1e–6) from genotype likelihoods using the SAMtools genotype likelihood model (Li et al. 2009). Allele frequencies were estimated assuming a fixed major but unknown minor allele (Kim et al. 2011) and a uniform prior. Sites with a minor allele frequency less than 0.05 and sites missing data in more than half of the individuals were excluded.

After establishing loci that met the above filtering criteria, for each individual, we sampled a single sequencing read at each locus on Omy5 (415 total), the chromosome with the migration-associated block of linkage disequilibrium (Pearse et al. 2014), and used the base call from that read as input data for downstream analyses (see below). In other words, for each individual, one read was sampled at each SNP passing the filtering step regardless of sequencing depth at that SNP. The single-read sampling approach mediates the effects of coverage differences between individuals and allows the inclusion of data from positions and samples with low coverage, facilitating the inclusion of a larger number of samples than would be possible with other approaches, such as calling genotypes (Korneliussen et al. 2014). The single-read base calls were used to generate an individual-by-loci matrix, where 1 represented the major allele and 0 represented the minor allele. Missing data were replaced with the mean allele call at the locus, which corresponds to the major allele frequency, to conduct DPAC, which requires a full numeric matrix (Jombart et al. 2011; see below).

To determine life history genotype groups, we conducted a DAPC (Jombart et al. 2010) with the matrix described above. Owing to its large size and high divergence, the variation in the Omy5 region dominates the discriminate analysis (Fig. 2). We used the “find.cluster” method implemented in R package “adegenet” (Jombart et al. 2011), a method that uses model selection to infer genetic groups by partitioning genetic variation into between- and within-group variation. We calculated Bayesian information criteria (BIC) for cluster models including k = 1 to k = 10 clusters and calculated the decrease in BIC between models to identify the optimal number of clusters (Jombart et al. 2010), methods akin to choosing the number of clusters in STRUCTURE (Evanno et al. 2005). Missing data were replaced with the mean value (major allele frequency) at each locus, which led to individuals with large amounts of missing data being grouped into the heterozygote group. To avoid false heterozygotes, we included individuals who had data at a minimum of 165 SNPs (assignment to the heterozygote group increases substantially among individuals who are missing more than this number of SNPs; see also Appendix A). This decision resulted in excluding life history genotype classifications for 54 individuals, 41 of which were assigned heterozygote-genotypes.

We conducted additional analyses on a subset of samples to determine their genotypic sex using presence–absence of a Y chromosome-linked marker, using sequences described by Brunelli et al. (2008). We used 1–10 ng of DNA as a template for a Taqman SNP genotyping assay (Life Technologies Corporation, Carlsbad, Caifornia). PCRs were done in a 10 μL volume containing Taqman GT Master Mix and custom Taqman probes and primers for OMY1–2SEXY and an autosomal locus to distinguish between a lack of template and lack of Y chromosome. The amplification was conducted on a QuantStudio3 (ThermoFisher Scientific) and consisted of a 10 min hold at 95 °C and then 40 cycles of 15 s at 95 °C and 1 min at 60 °C. Each plate included one male and one female control in addition to two blanks. Control samples for each sex were collected from known-sex adult O. mykiss at the Warm Springs Hatchery in Santa Rosa, California. We also calculated the sex ratios for a random subset of juvenile, age-0 fish (<85 mm FL, n = 239) to confirm a 1:1 sex ratio at the juvenile life stage.

We assigned observed life history ecotypes based on mark–recapture histories of individual fish and body size. In brief, individuals were assigned “migratory” if they were detected moving downstream during the spring smolt outmigration window and were assigned “resident” if they exceeded a size threshold (described below). To detect downstream movement, we installed stationary antennas at the mouths of each creek. We used multiplex readers from Oregon RFID (Portland, Oregon), and antennas were operated from November 2014 through May 2018 in Elder and Fox creeks (installation November 2014 in Elder and May 2015 in Fox). Antennas were located 200 and 350 m upstream of the mouth in Elder Creek and 175 and 195 m upstream of the mouth in Fox Creek. We attempted to operate antennas continuously, but high-flow events dislodged the antenna wires for periods during the wet season (up to 19 consecutive days). Here we focus on detections during the smolt migration window, from February to May, which is the period when O. mykiss smolts have been observed migrating in the Eel River watershed (Brown 1990). Individuals were assigned a “migratory” ecotype if they were detected moving downstream past the antenna arrays during these months (first detected at the upstream antenna and then at the downstream antenna). A subset of fish (n = 30) were removed due to detection histories that suggested persistent local movement and (or) shed tags (i.e., detected moving both upstream and downstream and were detected at the antenna over a span of time > 36 h, with 9–6571 detections per tag). There was a subset of individuals who were only detected once in the migration time frame, which only occurred when one antenna was functioning. Because all of the directional detections (n = 98) that occurred during the smolt outmigration period were in the downstream direction, we assumed that these single detections also represented fish moving downstream (n = 175). It is important to note here that we are assigning individuals to “migratory” ecotypes, rather than “anadromous”, and some individuals who are moving out of these headwater tributaries may be migrating to another part of the watershed rather than to the ocean.

We assigned individuals to a resident ecotype if they were >160 mm FL in July based on several lines of evidence. First, no individuals larger than 160 mm were detected moving downstream past the antennas during the outmigration window (98.2% of fish detected were less than 155 mm FL). Second, mark–recapture data (Kelson and Carlson 2019) and length frequency plots (refer to online Supplementary material, Fig. S1¹) indicate that a fish of 160 mm FL in the summer months is typically 2+ years old. The dominance of age-0+ and age-1+ fish in the out-migrants has been noted previously from O. mykiss smolt trapping in the South Fork Eel River (Brown 1990) and from a smolt trap operated in nearby Pudding Creek (Mendocino County, California; Ohms et al. 2014). Finally, this size cut-off is larger (i.e., more conservative) than the size threshold of 150 mm that was applied in nearby populations of O. mykiss to assign “resident” fish (streams in central California, Rundio et al. 2012).

To determine if sex ratios differed from 1:1 in juvenile fish, out-migrating fish, and resident fish, we used exact binomial tests, computed in R (R Core Team 2017). We also compared the genotype frequency of migratory ecotype fish and resident ecotype fish against the “baseline” genotype frequency of juvenile fish (less than 85 mm FL) using a χ² test in R.

Next, we modeled the relationship between life history genotype, genotypic sex, and life history ecotype using generalized linear models in R. Specifically, we used binomial models, with migratory fish assigned a value of 1 and resident fish assigned a value of 0. We calculated BIC for six models with different combinations of predictor variables: (1) sex; (2) life history genotypes; (3) sex and life history genotypes; (4) sex, life history genotype, and their interaction; (5) sex, life history genotype, and sample location; and (6) sex, life history genotype, their interaction, and sample location. We also calculated r² using the “rsq” package (Zhang 2018) for models including sex, genotype, and the combination of the two (models 1–3 listed above) to compare the ability of these factors to predict life history expression independently and in combination.

We combine data from Fox and Elder creeks in the results below for two reasons. First, we had limited sample sizes from Fox Creek (n = 38 individuals out of a total of 284 assigned a life history ecotype were from Fox). Second, the best-fit generalized linear binomial model did not include capture location (see Results).

We used an equation from Ohms et al. (2014) to estimate the proportion of the combined population that migrates, given the sex ratios in the out-migrating fish, the resident fish, and the premigration aged fish. In this model, if sex ratios are skewed for only the resident fish, but close to 1:1 for the out-migrating fish, then the number of individuals expressing residency is a small proportion of the total population (see also Ohms et al. 2014). We applied the same equation, substituting migratory allele frequencies for sex ratios:

where P is the proportion of migrants, m is the migratory allele frequency of migratory ecotype fish, b is the baseline migratory allele frequency of premigration-aged fish, and r is the migratory allele frequency of resident fish.

When grouping individuals into clusters based on Omy5 SNPs, we found the largest reduction in BIC occurred between models with k = 2 and k = 3 clusters (Fig. S2¹), which aligned with our prediction of three groups (resident, heterozygous, and migratory) corresponding to each genotype. We used group membership (k = 3) to call individuals as resident, heterozygous, or migratory genotypes (Fig. 2). Herein, we refer to the group assignment as “genotype”, where the “migratory genotype” means homozygous for the haplotype associated with the migratory life history, “resident genotype” means homozygous for the haplotype associated with the resident life history, and “heterozygous” means one copy of each haplotype.

We obtained genotype data from n = 4303 out of 4332 (99.3%) fish. We then assigned life history genotypes to n = 3450 fish (80.1% of those that were genotyped), with individuals being excluded due to missing data at many SNPs on Omy5 (see above; Appendix A). See Table S1¹ for a summary table of life history genotypes by year and age class. Across 4 years of sampling (2014–2017), the overall genotype frequency for juvenile fish (<85 mm FL, n = 2485) in these streams was 29.6% migratory, 40.2% heterozygotes, and 30.2% resident. We successfully assigned life history ecotypes and genotypes to 284 fish, which are summarized by genotype and sex (Table 1), and the data set is available on Dryad Digital Repository (Kelson et al. 2019).

We found that migratory genotype fish were unlikely to remain in the streams as resident fish, as resident ecotypes were composed of 57.1% resident genotypes, 36.5% heterozygous genotypes, and only 6.3% migratory genotypes (Fig. 3). This frequency of genotypes differed from the frequency of genotypes in juvenile fish in a χ² test for given probabilities (χ² = 27.0, p < 0.001). Similarly, we found that the proportion of migratory alleles decreased in a cohort of fish from age-0 to age-2+ (i.e., as the cohort aged; Fig. 4a), also indicating that migratory genotype fish were less likely to remain in the stream as older fish.

In contrast, we found that resident genotype fish often expressed the migratory ecotype. The genotype frequencies in the migratory ecotypes (24.7% migratory, 45.3% heterozygote, and 30.0% resident genotypes) were very similar to the baseline juvenile genotype frequencies and did not differ in a χ² test (χ² = 2.5, p > 0.1; Fig. 3). Similarly, we found that migratory ecotype fish did not always include a higher proportion of migratory alleles when compared with all the fish caught the previous summer (the baseline for that year; Fig. 4b).

For juvenile fish (i.e., age-0, <85 mm FL, n = 252), 48.8% were male and 51.2% were female, which was not significantly different from a 1:1 ratio (binomial test, p = 0.75). For fish assigned a resident ecotype (n = 58 with sex data), 77.6% were male (binomial test, p < 0.01, 95% confidence interval from 64.7% to 87.5%). In contrast, migratory ecotype fish (n = 153 with sex data) were 38.1% male, which differed significantly from an even sex ratio (binomial test, p < 0.01, 95% confidence interval from 30.4% to 46.2%).

Migratory ecotypes consisted of fish from every life history genotype (migratory, resident, and heterozygous) but were female-biased (Fig. 3). In contrast, resident ecotypes consisted primarily of resident genotypes and were male-biased (Fig. 3). Notably, there were no female migratory genotype fish who expressed the resident ecotype, and correspondingly, there were very few male resident genotype fish who expressed the migratory ecotype (Table 1; Fig. 3).

We found that the best model (lowest BIC) describing whether individuals expressed a migratory or resident ecotype included life history genotype and sex (Table 2). The probability of out-migrating increased with the addition of migratory alleles; migratory genotypes were the most likely to out-migrate (probability for males = 0.91 and for female = 0.99), followed by heterozygotes (male = 0.62, female = 0.92), and then resident genotypes (male = 0.31, female = 0.77; Fig. 5a). Based on effect sizes of parameters, sex was more important than genotype in explaining variation in ecotype (Fig. 5b, with heterozygote–female as the intercept; z = 6.0, p < 0.01; migratory genotype: z = 2.4, p < 0.05; resident genotype: z = –3.2, p < 0.01; and male: z = –4.9, p < 0.01). Similarly, the correlation between genotype and life history ecotype (r² = 0.20) was lower than the correlation between sex and life history ecotype (r² =0.31), but including both factors had the highest correlation (r² = 0.45).

Using the sex ratio data from our system as input data to the equation from Ohms et al. (2014), we estimated that 69.9% of the population migrated. Using migratory allele frequencies in the migratory versus resident ecotypes, we estimated that 88.3% of the population migrated.

In our system, genetic sex was correlated with observed life history, with resident ecotypes dominated by males and migratory ecotypes dominated by females (Fig. 3, Fig. 5a). The male-biased sex ratio in large, resident fish that we observed (78%) was close to the ratio observed (83%) in Big Creek (Monterey County, California; Rundio et al. 2012). In contrast, approximately 1280 km north, in the South Fork John Day River (Oregon), there was a 1:1 sex ratio among older fish (Ohms et al. 2014). This variation among systems suggests that the propensity for male bias in resident fish may depend on the environmental context. Overall, favorable freshwater growth conditions lead to a higher likelihood of freshwater maturation in male O. mykiss (McMillan et al. 2007; Doctor et al. 2014; Kendall et al. 2014), with males tending to be more plastic than females in exploiting freshwater resources in partially migratory populations (Berejikian et al. 2014). The pattern of male maturation in high growth conditions has been observed across salmonid systems (Jonsson and Jonsson 1993; Hendry et al. 2003), including widespread precocial maturation of age-0 male S. salar following a flood-induced pulse of food resources (Letcher and Terrick 1998). Additionally, the decision to mature in fresh water for males may be driven by access to mates (Gross 1985; Fleming and Reynolds 2003; Sloat et al. 2014).

Similarly, female bias in migrating fish also has been documented in salmonid systems. For example, Ohms et al. (2014) reported female-bias is O. mykiss out-migrants, with females representing 56%–76% of out-migrants in streams from northern California, Washington, Oregon, and Idaho, comparable to our results (62% females in the out-migrants). Long-term studies on Atlantic salmon also reported that out-migrants can be female-biased; across 7 years in the Saint Marguerite River, Quebec, Canada, females represented 50%–64% of the out-migrants (overall mean = 59%; Páez et al. 2011). Similarly, females represented 51%–76% (overall mean = 64%) of out-migrating S. salar across 11 years in the River Imsa, Norway (Jonsson et al. 1998). This pattern is consistent with the theory that females are more likely to benefit from migration due to enhanced growth opportunities and, hence, fecundity (Jonsson and Jonsson 1993; Fleming and Reynolds 2003; Hendry et al. 2003). Together these studies suggest that across salmonids, females are more likely to undertake ocean migration than males, but, similar to male bias maturation, this tendency may vary among populations (experiencing different growth conditions) or within populations (experiencing temporal variation in growth conditions).

It is important to note that male bias in resident fish does not necessarily imply female bias in out-migrants. When one life history strategy makes up a relatively small proportion of the overall population, sex ratios of the two life history forms become decoupled (Ohms et al. 2014). In our system, for example, the female skew we observed in out-migrants (62% female) was not as extreme as the male skew observed in freshwater residents (78% male). That sex ratios can be decoupled highlights the importance of estimating sex ratios in both the out-migrants and residents to understand the propensity of each sex to migrate within a given system and accordingly the utility of genetic sex for predicting life history ecotype.

We found that a model including both life history genotype and sex received the most support in terms of explaining individual life history decision; however, 55% of the variation remained unexplained. The decision to migrate is often considered a threshold trait, with individuals who reach a large enough body size by a certain time out-migrating the following spring (Satterthwaite et al. 2009, 2012). Beyond growth rates, the propensity to migrate can be influenced by body condition (Sloat and Reeves 2014), epigenetic regulation (Baerwald et al. 2015), and indirectly by maternal history (Liberoff et al. 2014), all of which we did not quantify in this study. Further research could identify if gene expression differs at the Omy5 loci for resident versus migratory ecotypes. Given that migration is a partially plastic trait, it is not surprising that life history strategy did not correlate perfectly with genotypic assignment in partially migratory populations.

Migratory ecotype individuals included a mixture of genotypes that was similar to the juvenile, or “baseline”, frequencies, including both resident and heterozygous genotype fish. This result emphasizes the importance of maintaining resident genotype fish in partially migratory populations because they are a source of migrants. Conservation of genetic diversity may buffer phenotypic variability in partially migratory populations. This result builds on a suite of earlier studies using otolith microchemistry to reveal that both resident and anadromous mothers can produce anadromous offspring (Zimmerman and Reeves 2000; Riva-Rossi et al. 2007; Zimmerman et al. 2009; Hodge et al. 2016). Similarly, studies from a population of O. mykiss isolated above a barrier for over 70 years revealed that the fish are still able to undergo smoltification, although with reduced fitness (Thrower and Joyce 2004; Thrower et al. 2004). In general, this body of work emphasizes the plastic nature of life history expression in partially migratory populations of O. mykiss, further demonstrating the importance of conserving the full suite of diversity.

A vexing management problem for partially migratory populations such as O. mykiss is determining the proportion of the population expressing the migratory ecotype of conservation concern. In this study we present two genetic tools, genotyping at Omy5 and sex-genotyping, both of which are useful for understanding migration behavior within partially migratory populations. We suggest that genotyping at Omy5 is more useful at the population level (as demonstrated by Abadía-Cardoso et al. 2011; Pearse et al. 2014; Apgar et al. 2017) than at the individual level given that 55% of the variation in life history expression was unexplained by genetic sex and life history genotype. Our results indicate that genetic sex provides more explanatory power than life history genotype at the individual level. Furthermore, determining genetic sex is a relatively easy and inexpensive procedure and can be done using nonlethal tissue samples and qPCR to genotype at two SNPs (one autosomal control and one Y-chromosome SNP).

Individual level life history decisions scale up to population level patterns, and quantifying these patterns may be of interest to managers seeking to restore or monitor anadromy in partially migratory populations. Using genetic tools, we generated two estimates for the percentage of the population expressing migration, the first based on sex ratios (estimated 70% migration) and the second based on migratory allele frequencies (estimated 88% migration) in these streams. Both approaches suggest that most of the population migrates from these streams. The mathematical models used to generate these estimates included simplifying assumptions such as that of equal mortality rates between sexes or genotypes before the life history decision window or before the samples are collected. We suggest that using the sex ratios, rather than migratory allele frequency, may be less likely to violate these assumptions. This is because other traits that have been linked to Omy5, including development rate (Miller et al. 2012) and growth (Nichols et al. 2008), could be related to age-specific mortality. Future studies could test the assumptions of these models by estimating sex-specific and genotype-specific juvenile mortality rates in fresh water for wild O. mykiss to further evaluate the utility of this approach for estimating the proportion of the population that is migratory. If these assumptions are robust, then determining the sex ratios of different groups in a population (smolts versus residents) may be a cost-effective way to estimate the proportion of a population that out-migrates. Furthermore, this tool may be useful for understanding population-level responses in migration behavior to restoration projects (e.g., dam removals, habitat improvements projects).

Within partially migratory populations of O. mykiss and other salmonid fishes, there is interest in conserving and restoring migration where it has been lost. Evaluating the success of conservation and restoration efforts requires an understanding of the proportion of the population expressing migration, which can be tackled at the population level (by estimating genotype frequencies or the proportion of fish expressing each life history) or the individual level (by correlating life history ecotype with life history genotype, as we did here). Our results suggest that genetic tools, especially when combining information on genetic sex and life history genotype, can be useful in estimating life history ecotype in partially migratory populations. However, 55% of the variation in life history ecotype remained unexplained after accounting for genetic sex and life history genotype. Our results suggest that using these tools at the population level is more appropriate for management decisions because life history genotype does not perfectly predict life history expression at the individual level. More generally, our results emphasize the importance of conserving genetic diversity in protected, partially migratory populations, including conserving resident fish in addition to migratory fish, because resident genotypes can give rise to migratory ecotypes.