Uncertainty about the magnitude, frequency, location, and timing of the nonlocal harvest of sockeye and chum salmon was the impetus for the Western Alaska Salmon Stock Identification Program. The program was designed to use genetic data in mixed stock analysis to reduce this uncertainty. A baseline of allele frequencies in spawning populations is required for use in mixed-stock analysis to estimate the stock of origin of harvested fish. This report describes the methodology used to understand the population genetic structure among chum salmon populations and to build and test a baseline for use in mixed stock analysis of chum salmon. Of the 35,921 fish from 434 collections selected to be genotyped, the final baseline was composed of 32,817 fish from 402 collections representing 310 populations. Average population sample size was 106 fish. Reporting groups were determined through a combination of stakeholder needs and identifiability using genetic information, as measured using proof tests. The final reporting groups included Asia, Kotzebue Sound, Coastal Western Alaska, Upper Yukon River, Northern District (Alaska Peninsula), Northwest District (Alaska Peninsula), South Peninsula (Alaska Peninsula), Chignik/Kodiak, and East of Kodiak.

Tissue Sampling

Baseline collections

Baseline samples were collected from spawning populations of chum salmon ranging from Korea to the State of Washington (Table 2). Tissue samples were collected by ADF&G and collaborators. In some cases, only DNA extracts were available, especially for collections from Asia. Target sample size for baseline populations was 95 individuals to achieve acceptable precision for estimating allele frequencies (Allendorf and Phelps 1981; Waples 1990a) and to accommodate our SNP genotyping platform. It was evident early on that a sample size of 95 fish was not possible for some collections.

Selection of baseline collections to genotype

We selected a subset of collections to include in the WASSIP baseline to represent 1) population abundance, 2) genetic diversity, 3) geographic coverage, and 4) among-year variation of allele frequencies within populations. We attempted to find the greatest representative value by balancing maximum representation and the high cost of choosing fish from every collection with minimal representation and the low cost of choosing fish from a small subset of collections. Unlike the baseline for sockeye salmon, we chose to include populations from throughout the species' range, because past studies indicated that chum salmon from spawning areas throughout the Pacific Rim are caught in WASSIP-area fisheries (Seeb and Crane 1999a, 1999b; Seeb et al. 2004).

Laboratory Analysis

Developing and ascertaining SNPs for WASSIP

At the start of this project, ADF&G used a panel of 53 SNP markers for most chum salmon analyses. Because of the need for greater resolution among Western Alaskan stocks, we contracted the discovery of additional SNP markers for chum salmon with the International Program for Salmon Ecological Genetics (IPSEG; http://fish.washington.edu/research/seeblab ) at the University of Washington. These efforts were based on cDNA sequences from two chum salmon from the Susitna and Delta rivers (Seeb et al. 2011a). This SNP development added 37 SNPs to those already available for chum salmon for use in WASSIP. Subsequent rounds of SNP development at the University of Washington were based on 16 fish from 4 populations from Western Alaska and increased the total number of SNPs to 198 (Petrou 2012). TaqMan™ assays were developed, or were available, for all 198 SNPs, including the original 53-SNP set.

We used a combination of assessments (Figure 1) to select the 96 most useful SNP markers for WASSIP (see below). These measures were performed on genotypic data collected for 188 SNP markers on 30 populations. This dataset was chosen to include populations relevant to WASSIP, and an emphasis was placed on choosing populations from Western Alaska. We refer to this dataset as the backbone marker set. This process of marker selection is described in DeCovich et al. (2012b), but a brief description is given here.

We eliminated markers from further consideration if they did not conform to Hardy-Weinberg Expectations (HWE). Markers which showed significant genotypic association with each other (linkage disequilibrium) were combined into composite genotypes and tested for significant increase in information when treated together (DeCovich et al. 2012b). If there was no significant increase in information, the least useful of the associated markers was dropped from further consideration. If there was a significant increase in information the combined marker was retained, referred to hereafter as a composite locus and SNP markers that are not associated with any other marker are referred to as loci.

Scored assessments were used to determine the relative usefulness of backbone loci for MSA performance with weights set primarily to select useful markers for two objectives: 1) to increase MSA performance to distinguish among Western Alaska regions, and 2) to maintain MSA performance to distinguish among regions throughout the rest of the species' range (Figure 2). Weights were assigned to measures based on perceived benefit of the measure to achieve the objectives and the highest weights were placed on measures associated with the first objective (Figure 3). For each measure, j (j = 1, 2, ..., J), raw values for each backbone locus, l (l = 1, 2, ..., L), were linearly scaled into scores between 0.0 (lowest) and 1.0 (highest) using the equation:

S'_l,j = (S_l,j - S_min,J) / (S_max,j - S_min,j)

For measure j, S_l,j is the raw value at backbone locus l, S_min,j is the minimal value across all backbone loci, S_max,j is the maximal value across all backbone loci, and S'_l,j is the scaled score for backbone locus l. Scores for each locus were summarized across measures using the following equation:

S_l = ∑_j=1-J (w_j S'_l,j) / (∑_l=1-L S'_l,j)

Here w_j is the weight assigned to measure j , and S_l is the final score given to backbone locus l. These scores were then ranked from highest to lowest. The 96 highest-ranking backbone loci were then assessed for laboratory performance. If a marker performed poorly in the laboratory, it was replaced by the next-highest-ranking backbone locus.  The process was repeated for each poorly performing backbone locus, until a final set of 96 backbone loci were selected.

Assaying genotypes

We extracted genomic DNA from tissue samples using a DNeasy^® 96 Tissue Kit by QIAGEN^® (Valencia, CA). A multiplexed preamplification PCR of 96 screened SNP markers was used to increase the concentration of template DNA. Reactions were conducted in 10 μL volumes consisting of 4 uL of genomic DNA, 5 μL of 2×Multiplex PCR Master Mix (QIAGEN) and 1 μL each (2 μM SNP unlabeled forward and reverse primers). Thermal cycling was performed on a Dual 384-Well GeneAmp^® PCR system 9700 (Applied Biosystems) at 95°C for 15 min followed by 20 cycles of 95°C for 15 s, 60°C for 4 min, and a final extension at 4°C. The preamplified DNA was then loaded into a Fluidigm^® 96.96 Dynamic Array in a post-PCR laboratory at ADF&G. The Fluidigm^® 96.96 Dynamic Array contains a matrix of integrated channels and valves housed in an input frame. On one side of the frame are 96 inlets to accept the sample DNA from individual fish and on the other are 96 inlets to accept the assays for 96 SNP markers. Once in the wells, the components are pressurized into the chip using the IFC Controller HX (Fluidigm). 

The 96 samples and 96 assays are then systematically combined into 9,216 parallel reactions. Each reaction is a mixture of 4 µl of assay mix (1×DA Assay Loading Buffer (Fluidigm), 10×TaqMan^® SNP Genotyping Assay (Applied Biosystems), and 2.5×ROX (Invitrogen)) and 5 µl of sample mix (1×TaqMan^® Universal Buffer (Applied Biosystems), 0.05×AmpliTaq^® Gold DNA Polymerase (Applied Biosystems), 1×GT Sample Loading Reagent (Fluidigm) and 60–400 ng/µl DNA) combined in a 7.2 nL chamber. Thermal cycling was performed with an Eppendorf IFC Thermal Cycler as follows: 70°C for 30 min for “Hot-Mix” step, initial denaturation of 10 min at 96°C followed by 40 cycles of 96° for 15 s and 60° for 1 min. The Dynamic Arrays were read on a Fluidigm^® EP1^TM System or BioMark^TM System after amplification and scored using Fluidigm^® SNP Genotyping Analysis software.

Assays failing to amplify with the Fluidigm system were reanalyzed on the Applied Biosystems platform using the preamplified DNA. Each reaction on this platform was performed in 384-well reaction plates in a 5 µL volume consisting of 5-40 ng/μl of template DNA, 1×TaqMan^® Universal PCR Master Mix (Applied Biosystems), and 1×TaqMan^® SNP Genotyping Assay (Applied Biosystems). Thermal cycling was performed on a Dual 384-Well GeneAmp^® PCR System 9700 (Applied Biosystems) as follows: an initial denaturation of 10 min at 95°C followed by 50 cycles of 92°C for 1 s and annealing/extension temperature for 1 min. The plates were scanned on an Applied Biosystems Prism 7900HT Sequence Detection System after amplification and scored using Applied Biosystems’ Sequence Detection Software version 2.2.  

Genotypes produced on both platforms were imported and archived in the Gene Conservation Laboratory Oracle database, LOKI.

Laboratory quality control

We conducted a quality control analysis (QC) to identify laboratory errors and to measure the background discrepancy rate of our genotyping process. The QC analysis was a collaborative effort between ADF&G Gene Conservation Laboratory and the University of Washington IPSEG lab. The QC analyses were performed by staff not involved in the original genotyping. We applied 3 methods to the QC, depending on the type of collection and when it was genotyped. We termed these the Original, Database, and New QC methods. All QC DNA underwent the preamplification process as described above, unless otherwise noted.

The Original QC method was the method used for QC prior to beginning WASSIP. This method consisted of regenotyping 8% of the fish genotyped in the original project using the same DNA extraction for the same SNPs assayed in the original project. Discrepancy rates were calculated as the number of conflicting genotypes, divided by the total number of genotypes examined. These discrepancy rates describe the difference between original project data and QC data and are capable of identifying assay plate errors but cannot detect DNA extraction plate errors (rotations, etc.), since they were based on the same extractions. This QC method was implemented only when all tissue was exhausted during the original extraction or when the DNA was extracted off-site.

The Database QC method compared new and old genotypes for the SNPs common to our current and previous baselines (Jasper et al. 2012). After selecting the 96 loci to analyze chum salmon (DeCovich et al. 2012b), we assayed each collection for all 96 SNPs on a single chip. Since some of these SNPs had been assayed earlier in some of the collections, we were able to compare the apparent genotypes for these SNPs in these collections. Discrepancy rates were calculated as the difference between old and new genotypes for these SNPs. These comparisons identified genotyping errors, but could not detect DNA extraction errors, since they were based on the same extractions. When a difference appeared, we used the new genotype. This provided consistency between all of the data used in our analyses. Collections requiring either the New or Database QC methods were run with the same assay plate so that the New QC method could identify errors in the assay plate that the Database method could not.

The New QC method is our current QC method and consists of re-extracting DNA from 8% of project fish and genotyping them for the original SNPs. These discrepancy rates measure differences between the original and QC genotypes and are capable of identifying extraction, assay plate, and genotyping errors. This QC method represents the error rate in our current genotyping procedure.

For all QC methods, error rates in the original genotyping can be estimated as half the rate of discrepancy by assuming that the discrepancies among analyses were due equally to errors during the original genotyping and to errors during quality control, and by assuming that at least one of these assays produced the correct genotype.

Statistical Analysis

Data retrieval and quality control

We retrieved genotypes from LOKI and imported them into R (R Development Core Team 2010). Subsequent analyses were performed with R routines unless otherwise noted.  We performed 2 analyses prior to statistical analyses to confirm the quality of the data.

First, we removed from further analyses genotypes of fish that were missing a substantial number of genotypes for individual SNPs. We used an 80% rule (Dann et al. 2009) to exclude fish missing genotypes for 20% or more loci, because these individuals likely had poor-quality DNA. The inclusion of individuals with poor-quality DNA might introduce genotyping errors into the baseline and reduce the accuracy of MSA.

Second, we identified individuals with identical genotypes and removed them from further analyses. Identical genotypes can result from sampling or extracting the same fish twice. Identical was defined as pairs of fish with the same genotypes in 83% of SNPs (80 of the 96 backbone loci in this study). One individual from each duplicate pair was removed from further analyses.

Hardy-Weinberg expectations

After calculating allelic frequencies for each backbone locus, we tested observed genotype frequencies for each baseline collection for conformance to HWE at each marker by Monte Carlo simulation with 10,000 iterations using the adegenet package (Jombart 2008). We combined probabilities for each collection across backbone loci using Fisher’s method (Sokal and Rohlf 1995) and removed collections that departed significantly from HWE after correcting for multiple tests with Bonferroni’s method (α = 0.05) from subsequent analyses.

Pooling collections into populations

When appropriate, we pooled collections following a step-wise protocol to obtain more accurate estimates of population allele frequencies. First, we pooled collections from the same geographic location, sampled at similar calendar dates but in different years, as suggested by Waples (1990b). We used Fisher’s exact test of allele frequency homogeneity (Sokal and Rohlf 1995) and based our decisions on a summary across loci using Fisher’s method. Collections that were not significantly different (P > 0.01) were pooled. After pooling, we treated these final collections as populations in our analyses. Finally, we tested genotypic frequencies for conformance to HWE following the same protocol described above to ensure that our pooling was appropriate and that tests for linkage disequilibrium would not result in falsely positive results due to departure from HWE.

Process for defining reporting groups

Defining reporting groups for MSA was an iterative process that took into account the following criteria: 1) sociological needs (suggested by stakeholders and fishery managers), 2) genetic population structure (MSA potential), 3) adequate representation in the baseline (number of individuals and representative value of genetic variation within groups), and 4) the expected number of fish from a reporting group potentially within a mixture (Habicht et al. 2012). To evaluate potential reporting groups, we used the following metrics of these 4 factors as guidelines: 1) utility of information for fishery managers and stakeholders, 2) 90% correct allocation in tests of the baseline’s ability to allocate to reporting groups, 3) 400 individuals from enough collections to adequately represent the genetic diversity in a reporting group, and 4) an expected contribution to a given mixture of at least 5% (or 20 fish) for the 400-fish mixtures proposed for WASSIP. The definition of reporting groups heavily depended on information from the Testing reporting groups for MSA and identifying biases section described below.

Removal of collections from the baseline

We removed collections from further analysis if they met one of the following criteria: 1) a collection did not meet our minimum desired sample size of 40 individuals after pooling, or 2) poorly documented collections that did not pool with other well-documented collections putatively sampled at the same location and time of year. In these cases, the less-documented collections were excluded.

Linkage disequilibrium

We re-evaluated the data for linkage disequilibrium between each pair of nuclear backbone loci in each population to ensure that subsequent baseline and mixed stock analyses would be based on independent markers. We used the program Genepop version 4.0.11 (Rousset 2008) with 100 batches of 5,000 iterations for these tests. We summarized the frequency of significant linkage disequilibrium between pairs of SNPs (P < 0.05) and further investigated pairs that exhibited linkage in more than half of the populations.

As done previously during locus selection (DeCovich et al. 2012b), we either removed one of the linked backbone loci or combined the pair into a composite locus for further analyses, if the pattern of linkage continued to provide information useful for MSA. We used the optimal rate of correct assignment (f_ORCA; Rosenberg 2005) as our measure of information. f_ORCA assesses the rate of correct allocation of simulated individuals to defined reporting groups based upon the markers in question. Composite loci generally have higher f_ORCA values than the single markers that form them, because combinations of alleles for 2 or more composite loci can exist in more forms (9 possible) than 2 single loci (4 alleles for a pair of loci). Simple comparisons of f_ORCA values would always suggest combining linked pairs into composite loci. However, there is a cost associated with composite loci as estimates of frequencies for 8 phenotypically scored alleles are less precise than estimates of 1 allele frequency at each of 2 loci for a given sample size.

To account for this cost, and to ensure that we combined only SNP pairs that provided significantly more information than the single SNPs in question, we compared the difference between f_ORCA values of the composite locus and the single locus with the greater f_ORCA value in the pair (Δ = f_ORCA-pair - max(f_ORCA-single1 ,  f_ORCA-single2)). This difference (Δ) was our test statistic, but since we did not know its underlying distribution, we conducted a sampled randomization test (Sokal and Rohlf 2005).  We randomly selected 1,000 SNP pairs, calculated Δ for each pair to empirically define the test-statistic distribution, and set the 90^th quantile of the distribution as a critical value (Δ₉₀). We then either combined linked backbone loci into a composite locus, if Δ was greater than this critical value, or we dropped the locus with the lower f_ORCA value, if Δ was less than the critical value.

This investigation of linkage disequilibrium is identical to that described in DeCovich et al. (2012b). That analysis produced the final set of 96 backbone loci used in this analysis, but it was repeated here to check for differences in associations due to the increase in populations in the full baseline.

We use the term locus in place of the term backbone locus from here on to describe the final set of loci used for analyses.

Treatment of mtDNA

The decision was made a priori to retain the 3 mtDNA SNP loci during the locus-selection process. These 3 markers were combined into single composite phenotypic locus, Oke_Cr30_Cr386_ND3-69, for the baseline. These mtDNA loci are useful for the identification of Japan-origin chum salmon, and our treatment of these loci is identical to that of Seeb et al. (2011b).

Analysis of genetic structure

Analysis of temporal variance

We examined the among-year temporal variation of allele frequencies with a hierarchical, 3-level Analysis of Variance (ANOVA). This analysis was performed on all populations with multiyear collections, but excluded all collections that were removed following the collection removal outlined above under Pooling collections into populations. We treated temporal samples as subpopulations, as described in Weir (1996). This method allowed the quantification of the sources of total allelic variation and permitted the calculation of the between-collection component of variance and the assessment of its magnitude relative to the between-population component of variance. This analysis was conducted using the software package GDA (Lewis and Zaykin 2001).

Visualization of genetic distances

We visualized pairwise F_ST estimates among collections from the final set of independent loci estimated with the package hierfstat (Goudet 2006). We constructed consensus Neighbor-Joining trees from 1,000 bootstrapped trees by resampling loci with replacement to assess the significance of nodes in the tree across loci. We plotted the consensus tree with the FigTree program (Rambaut 2007). These trees provided insight into the genetic structure among populations.

Determination of reporting groups

At the March 2011 AP meeting, the AP recommended sets of reporting groups deemed useful for answering questions of stock composition for the fisheries sampled in the study. The groups were arranged hierarchically, with the finest scale (greatest number of groups) desired given first, then with alternative broad-scale groups. The fine-scale groups included: 1) Asia, 2) Kotzebue Sound, 3) Norton Bay/Shaktoolik/Unalakleet, 4) Nome/Port Clarence/Golovin/Elim, 5) Yukon River, 6) Kuskokwim River 7) Togiak/Nushagak, 8) Eastern Bristol Bay, 9) Northern District Peninsula, 10) Northwest District Peninsula, 11) South Peninsula, 12) Chignik, and 13) East of WASSIP. The broadest scale included 6 groups: 1) Asia, 2) Kotzebue Sound, 3) Arctic, Yukon, Kuskokwim (AYK), 4) Bristol Bay, 5) Alaska Peninsula/Chignik, and 6) East of WASSIP. As a starting point for testing the performance of these groups, we chose a set of 10 groups between the fine-scale and broad-scale groups, and used an iterative approach to either collapse or expand the 10 groups based on the results of the previous round of analysis (Table 3). We presented the results of these tests at the September 2011 Advisory Panel meeting.

Baseline evaluation for mixed stock analysis

We assessed the identifiability of regional reporting groups in mixtures. These tests were used to determine if the underlying genetic structure supported particular reporting groups for MSA. These tests also provided insights into potential allocation biases. The results of these tests provided key insights in interpreting the results from MSA of WASSIP mixtures.

To assess the identifiability of regional reporting groups in mixtures we conducted 100% proof tests, in which we sampled 200 individuals without replacement from each reporting group and analyzed them as a mixture against the reduced baseline. These tests provided an indication of the power of the baseline to produce accurate MSA estimates under the assumption that all the populations from a reporting group were represented in the baseline. The AP and TC set a guideline that correct allocation for these single-reporting group tests should exceed 90% to be considered adequate (Seeb et al. 2000).

BAYES protocol

Stock compositions of these test mixtures were estimated with the program BAYES (Pella and Masuda 2001). The Bayesian model implemented in BAYES uses a Dirichlet distribution with specified parameters as the prior distribution for the stock proportions. We defined prior parameters for each reporting group to be equal (i.e., a flat prior) with the prior for each reporting group subsequently divided equally to populations within that reporting group.  We set the sum of all prior parameters to 1 (prior weight), which is equivalent to adding 1 fish to each mixture (Pella and Masuda 2001). We ran 5 independent Markov Chain Monte Carlo (MCMC) chains of 40,000 iterations with different starting values and discarded the first 20,000 iterations to remove the influences of the initial start values. We defined the starting values for the first chain such that the first 1/5 of the baseline populations summed to 0.9 and the remaining populations summed to 0.1. Each chain had a different combination of 1/5 of baseline populations summing to 0.9. 

We combined the second halves of these chains to form the posterior distribution and tabulated mean estimates and 90% credibility intervals from a total of 100,000 iterations.  We also assessed the within- and among-chain convergence of these estimates using the Raftery-Lewis and Gelman-Rubin diagnostics, respectively. These values compare variation of estimates within a chain (Raftery and Lewis 1996) and the total variation among chains (Gelman and Rubin 1992), respectively. If the Gelman-Rubin diagnostic for any stock group estimate was greater than 1.2 and the Raftery-Lewis diagnostic suggested that each chain had not converged to stable estimates, we reanalyzed the mixture with 80,000-iteration chains following the same protocol. We repeated this procedure for each reporting group mixture. A critical level of 90% correct allocation was used to determine if the reporting group was acceptably identifiable (Seeb et al. 2000). We visualized these results as barplots using the gplots package (Warnes 2010).

Genepop and rubias.