DOI: 10.5061/dryad.hx3ffbgc2
title: The Pacific lamprey genomic divergence, association mapping, temporal Willamette Falls, spatial rangewide datasets

description: High rates of dispersal can breakdown coadapted gene complexes. However,
 concentrated genomic architecture (i.e., genomic islands of divergence)
 can suppress recombination to allow evolution of local adaptations despite
 high gene flow.  Pacific lamprey (Entosphenus tridentatus) is a
 highly dispersive anadromous fish.  Observed trait diversity and
 evidence for genetic basis of traits suggests it may be locally adapted.
 We addressed whether concentrated genomic architecture could influence
 local adaptation for Pacific lamprey. Using two new whole genome
 assemblies and genotypes from 7,716 single nucleotide polymorphism (SNP)
 loci in 518 individuals from across the species range, we identified four
 genomic islands of divergence (on chromosomes 01, 02, 04, and 22). We
 determined robust phenotype-by-genotype relationships by testing multiple
 traits across geographic sites. These trait associations likely explain
 genomic divergence across the species’ range. We genotyped a subset of 302
 broadly distributed SNPs in 2,145 individuals for association testing for
 adult body size, sexual maturity, migration distance and timing, adult
 swimming ability, and larval growth.  Body size traits were
 strongly associated with SNPs on chromosomes 02 and 04. Moderate
 associations also implicated SNPs on chromosome 01 as being associated
 with variation in female maturity. Finally, we used candidate SNPs to
 extrapolate a heterogeneous spatiotemporal distribution of these predicted
 phenotypes based on independent datasets of larval and adult collections.
 These maturity and body size results guide future elucidation of factors
 driving regional optimization of these traits for fitness. Pacific lamprey
 is culturally important and imperiled.  This research addresses
 biological uncertainties that challenge restoration efforts.
description: The population genomic
 divergence dataset of 7716 quality-filtered SNP loci. <br> Divergence
 mapping<br> For characterization of SNP densities and FST
 statistics, we used a set of 7,716 unique SNP loci from previously
 published RAD-seq datasets (Hess et al. 2013; Smith et al. 2018), which
 passed a set of population genetic quality control filters (Supplemental
 Materials).  This set of 7,716 unique SNPs was a combination of
 overlapping groups of SNPs from a previous dataset (Hess et al. 2013; SNPs
 N = 8,772) and a de novo linkage mapping dataset (Smith et al. 2018; SNPs
 N = 7,977).  BOWTIE2 (Langmead and Salzberg 2012) was used to align these
 two datasets to the male reference assembly to define homologous loci.
  For the 7,716 total SNPs passing the QC filters, 4,046 loci were unique
 to Hess et al. 2013, 1,418 loci were unique to Smith et al. 2018, and
 2,252 SNPs were shared across datasets.  Marker positions based on BOWTIE2
 alignments were compared between Pacific lamprey male and female genomes
 and the Pacific lamprey male and sea lamprey (Petromyzon marinus) male
 gametic genome (GenBank assembly accession: GCA_002833325.1) to
 characterize synteny.<br> Using these 7,716 SNPs genotyped for the
 same individuals from Hess et al. (2013; i.e., 16 collections with
 &gt;20 individuals which totaled 482 individuals; Table 1), LOSITAN
 (Antao et al. 2008) was run using parameter settings of 50,000
 simulations, confidence interval of 0.99, false discovery rate set to 0.1,
 subsample size of 20, simulated FST of 0.019 and an attempted FST of
 0.021. We considered loci candidates for positive selection (adaptive
 loci) above a probability level of 0.995, and neutral loci were defined as
 falling between the 10th and 90th quantiles of the FST distribution. Any
 remaining SNPs were conservatively considered undetermined (neither
 adaptive nor neutral).<br> ---- The six datasets used for
 association mapping with 302 unique SNP loci.
 Association testing<br> Genotyping-in-thousands by
 sequencing (GT-seq, Campbell et al. 2015) was employed to genotype 308
 genetic markers for the association testing analyses.  The GT-seq 308 loci
 were a subset of markers developed from the paired end consensus reads
 from the Hess et al. (2013) RAD-seq dataset.  The selection of loci and
 steps in development are described in detail in Supplemental Materials.
  Locus selection began with a group of 457 total SNP loci considered in
 round 1, which included 120 for which TaqMan assays had already been
 designed (Hess et al. 2015).  Final optimization left 308 loci that worked
 best in GT-seq genotyping.  For all datasets used below in the association
 testing we filtered out individuals missing &gt;10% of genotypes at
 the 308 loci.  Excluding the four species diagnostic loci and two
 duplicated loci provided 302 unique loci for association tests.<br>
 There were six datasets, five comprised of adults (JDD, S_BON, T_BON,
 WFA♀, and WFA♂) and one comprised of larvae (GAR), with which we performed
 association testing (Table 1).  Adult datasets were from the following
 three locations: males (WFA♂, N=136) and females (WFA♀, N=133) from
 Willamette Falls collected in 2016 (Willamette River, Oregon City, OR;
 205.6 Rkm upstream from the Columbia River mouth), two datasets (S_BON,
 N=295 and T_BON, N=883) from Bonneville Dam in 2014 (235.1 Rkm upstream
 from the Columbia River mouth), and one dataset (JDD, N=656) from John Day
 Dam in 2014 and 2015 (346.9 Rkm upstream from the Columbia River mouth).
  The following five adult traits were measured on all adult datasets:
 ordinal “day” of collection (timing of migration to the sample point),
 girth (mm), total “length” (mm), weight (g), and distance between dorsal
 fins (“interdorsal”, mm).  Interdorsal measurements have been suggested to
 serve as an indicator of maturation status in Pacific lamprey because the
 distance tends to decrease with maturation (Clemens et al. 2009).  We
 measured an additional migration trait for three adult datasets (S_BON,
 T_BON, and JDD) via a combination of passive integrated transponder (PIT)
 and radio tagging of individual fish and observing their furthest upstream
 detection from the release location (“Rkm”).  Further, since the males and
 females collected at Willamette Falls (WFA♂ and WFA♀) were being
 harvested, we were able to measure gonad weight as a proxy for maturity in
 those datasets.  Finally, a subset of the adult dataset from Bonneville
 Dam (S_BON) was used in a swim trial experiment within a flume (Kirk et
 al. 2016), in which the following three swimming behavioral traits were
 measured: “approached” experiment, passed challenge (“pass”), and passed
 challenge without fallback (“passrep”). Details of these swimming
 performance experiments can be found in Kirk et al. (2016) and
 Supplemental Materials.<br> A single group of larvae were
 artificially propagated using adults captured at Bonneville Dam.  These
 larvae were reared in a common garden experiment to generate early larval
 growth (“GAR”) rate data (N=337).  All larvae were spawned in the spring
 of 2015 and allowed to rear from 30 to 163 days after hatching at which
 point they were measured for growth. Growth rate was measured as length /
 time (“growth”), and also corrected growth rate [“growth rate_b”; (length
 – 4 mm) / time] to correct for length at hatch (~4 mm).<br>
 Intercorrelation among all measured traits in these six datasets (i.e.,
 JDD, S_BON, T_BON, WFA♀, WFA♂, and GAR) was examined (based on Pearson’s
 r) to avoid excessive redundancy of predictor variables (│r│ &gt;
 0.95), and P-values were calculated (SAS Institute, Inc. 2000).  We
 performed univariate analyses using a general linear model (GLM) and a
 mixed linear model (MLM) with TASSEL v. 5.1.0 (Bradbury et al. 2007). The
 GLM is a fixed effects linear model that is used in TASSEL to identify
 significant associations between phenotypes and genotypes.  TASSEL takes
 population structure into account by using genetic principal coordinate
 axes as covariates in the model. The MLM is similar to GLM but includes
 both fixed effects (e.g. population structure, and genetic marker) and
 random effects (i.e., relationships among individuals) and can thus
 account for both population structure and kinship to reduce false positive
 associations (Yu et al. 2006).  Details on the covariates and ways in
 which loci were used taking population structure and relatedness into
 account in the GLM and MLM tests are provided in the Supplemental
 Materials. To account for multiple tests, only those associations with
 P-values less than the critical value as determined using the false
 discovery rate procedure described by Benjamini and Hochberg (1995) were
 considered significant.  The Benjamini and Hochberg (1995) false discovery
 rate approach has more power to detect significant differences than
 sequential Bonferroni correction (Narum 2006).  Critical values were
 calculated using the function p.adjust within the R package stats (RDC
 Team 2019). Association mapping<br> The 308 SNP
 loci in the GT-seq panel were aligned to reference genomes using BOWTIE2.
 There were 306 that were each assigned to a single location on the Pacific
 lamprey male genome (99.4%), covering 70 different chromosomes with an
 average of 4.4 loci per chromosome (range 1 – 22).  Marker locations were
 based on the alignments of marker sequences to the Pacific lamprey male
 and female genomes, homologous scaffolds of the sea lamprey genome, and
 positions on the previously published Pacific lamprey linkage map (Smith
 et al. 2018).<br> Adjusted P-values from the association testing
 described above were log transformed (-LOG10) and plotted by consensus
 genome position on the Pacific lamprey male genome. We tested correlation
 of association tests -LOG10(P) with FST from the rangewide divergence to
 understand whether trait associations may explain the high divergence
 observed at the rangewide scale for the subset of markers shared between
 datasets.  Among the 308 SNPs, there were 230 neutral SNPs, 41 adaptive
 markers SNPs, and a set of 31 “intermediate” SNPs that did not fit
 definitions of putatively neutral and putatively adaptive (divergence
 mapping).  Finally, four loci were species-diagnostic, and 2 loci were
 duplicated.  Therefore, there were 302 unique markers available for these
 association analyses. These markers included 38 SNPs that were mostly
 adaptive loci that were categorized into the following 4 groups of
 statistically linked loci: A (N=10), B (N=13), C (N=7), and D (N=8, Hess
 et al. 2013).<br> ---<br> Two temporal datasets of Willamette
 Falls (2014-2015) adults with 302 unique SNP loci
 <br> For the temporal datasets, we used individuals
 collected from two successive spawning runs at Willamette Falls (2014 –
 2015; N of 868 and 581, respectively) over which it was possible to
 randomly sample the majority of the annual adult migration of Pacific
 lamprey (typically February – August) in weekly strata. A daily abundance
 estimate (Whitlock et al. 2019) was used to expand candidate SNP genotypic
 proportions in the weekly strata.  We estimated the abundance and 95%
 confidence intervals of the candidate genotypes using a bootstrapping
 method (Steinhorst et al. 2017) that was automated in R for a broad set of
 fisheries applications that require stratified sampling (Thomas Delomas,
 PSMFC/IDFG, https://github.com/delomast/fishCompTools)".  One
 biological complexity was that a portion of the adults encountered before
 May probably overwintered and experienced shrinkage in body size due to
 advanced maturation (Beamish 1980). Therefore, in addition to
 characterizing allele frequencies of candidate SNPs Etr_464 and Etr_5317,
 we categorized fish by body length to provide insight into the transition
 between overwintered fish and newly-arrived migrants.<br>
 ---<br> The spatial dataset of larvae/juveniles rangewide with 85
 SNP loci<br> For the spatial dataset, we primarily used collections
 of larvae and juveniles (95% of dataset of N=3,435) but included some
 adult collections that were distributed widely across the species’ range.
  Larvae and juveniles were the ideal life stage to represent genotypic
 distributions of individuals that successfully spawned at discrete
 locations throughout the range. Adult collections were used to fill in
 portions of the range where larval samples were not available.  Genotyping
 was partially conducted with a TaqMan assay panel (Hess et al. 2015),
 which overlapped the GT-seq panel by 85 SNPs they had in common.  COLONY
 v. 2.0.6.5 (Jones and Wang 2010) was used to reconstruct full-sibling
 families (Wang 2004) using the 85 shared SNPs on each of the 70
 collections. We analyzed all collections together as one using the
 following parameter settings: polygamous mating for males and females
 without inbreeding, full-likelihood, medium length of run, no allele
 updating, and no sibship priors. Only 1 collection out of the 7 adult
 collections had full siblings (N=13, Stamp River, B.C.) which were
 maintained to accurately represent this small spawning segment.  We
 excluded duplicate genotypes, 797 full siblings, and collections with
 fewer than 5 individuals, resulting in a final set of 57 collections
 consisting of a total of 2,581 individuals each representing a unique
 family (Table S1). This dataset was then used to calculate allele
 frequencies across collections for the representative candidate SNPs
 Etr_464 and Etr_5317 within the adaptive regions on chromosomes 01 and 02,
 respectively.<br>  
description: There are 4 separate data files contained within this DOI, and
 the notes below are separated into four sections to accompany these
 datafiles.<br> The population genomic divergence dataset of 7716
 quality-filtered SNP loci<br> This is an excel file that includes
 tabs with an input file in genepop format that was used to run the FST
 outlier analysis. In addition, the metadata include the approximate
 genomic positions of the loci used to produce Manhattan plots for
 divergence mapping.<br> ---<br> The six datasets used for
 association mapping with 302 unique SNP loci.<br> The excel
 spreadsheet includes a README file to guide users through the tabs
 contained in the spreadsheet which include all the data necessary to run
 the association tests and the genome alignment positions to perform the
 association mapping.<br> ---<br> Two temporal datasets of
 Willamette Falls (2014-2015) adults with 302 unique SNP loci<br>
 This Excel spreadsheet contains all the necessary data to replicate our
 analysis of temporal distribution of the genotypes in the Willamette Falls
 adult Pacific lamprey collected in 2014 and 2015.
 Spread sheet contents:<br> WFA14_15RawData: The genotypic
 data collected from 295 loci and the individual trait data (length,
 capture date)<br> WFA14_15InputData: Only the subset of data used
 for the analysis of the study, i.e., the length category + Etr_5317
 genotypes, the Etr_464 genotypes, the Etr_5317 genotypes, and the capture
 date and weekly strata.<br> R_script_CI_bootstrapping: The exact R
 scripts run to perform the abundance estimates and 95% Confidence
 Intervals using the input data and 10,000 bootstraps<br> .csv files:
 The exact content in each csv file that is called out in the R scripts;
 Each run year (2014 and 2015) was divided in half by approximately 50% of
 the total estimated daily abundance (Early vs Late) and each category was
 estimated within each half of the run:<br> 2 run years x 2 halves of
 the run X 3 categories = 12 analyses (12 lines of R script to run total
 analysis).<br> Results: All abundance estimates and 95% CI for the
 categories of the 12 analyses and the charts used for the Figure
 S7 <br> ---<br> The spatial dataset of larvae/juveniles
 rangewide with 85 SNP loci<br> This Excel spreadsheet contains all
 the necessary information to recreate the metadata used to generate the
 Table S1 that appeared in the Supplemental Materials.<br> There are
 separate tabs within the spreadsheet containing the input files to run
 COLONY can be found on tabs "MakerTypeErrorRate" and
 "OffspringGenotype" which can be run as text files in the
 program COLONY.