Conservation genomics of anadromous Atlantic salmon
across its North American range: outlier loci identify the
same patterns of population structure as neutral loci
JEAN-S
EBASTIEN MOORE,*
1
VINCENT BOURRET,*†
1
M
ELANIE DIONNE,‡ IAN BRADBURY,§
PATRICK O’ REILLY,¶ MATTHEW KENT,** G
ERALD CHAPUT†† and LOUIS BERNATCHEZ*
*Institut de Biologie Int
egrative et des Syst
emes, Universit
e Laval, 1030 Avenue de la M
edecine, Qu
ebec, Qu
ebec G1V 0A6,
Canada, †Direction de la Protection de la Faune, Minist
ere des For
^
ets, de la Faune et des Parcs, Qu
ebec, Qu
ebec G1S 4X4,
Canada, ‡Direction de la Faune Aquatique, Minist
ere des For
^
ets, de la Faune et des Parcs, Qu
ebec, Qu
ebec G1S 4X4, Canada,
§Science Branch, Fisheries and Oceans Canada, 80 East White Road, St. John’s, Newfoundland A1C 5X1, Canada, ¶Science
Branch, Fisheries and Oceans Canada, Bedford Institute of Oceanography,1 Challenger Drive, Dartmouth, Nova Scotia, B2Y
4A2, Canada, **Centre for Integrative Genetics (CIGENE), Department of Animal and Aquacultural Sciences (IHA), Norwegian
University of Life Sciences, PO Box 5003, 1432 Aas, Norway, ††Fisheries and Oceans Canada, PO Box 5030, Moncton, New
Brunswick E1C 9B6, Canada
Abstract
Anadromous Atlantic salmon (Salmo salar) is a species of major conservation and
management concern in North America, where population abundance has been declin-
ing over the past 30 years. Effective conservation actions require the delineation of
conservation units to appropriately reflect the spatial scale of intraspecific variation
and local adaptation. Towards this goal, we used the most comprehensive genetic and
genomic database for Atlantic salmon to date, covering the entire North American
range of the species. The database included microsatellite data from 9142 individuals
from 149 sampling locations and data from a medium-density SNP array providing
genotypes for > 3000 SNPs for 50 sampling locations. We used neutral and putatively
selected loci to integrate adaptive information in the definition of conservation units.
Bayesian clustering with the microsatellite data set and with neutral SNPs identified
regional groupings largely consistent with previously published regional assessments.
The use of outlier SNPs did not result in major differences in the regional groupings,
suggesting that neutral markers can reflect the geographic scale of local adaptation
despite not being under selection. We also performed assignment tests to compare
power obtained from microsatellites, neutral SNPs and outlier SNPs. Using SNP data
substantially improved power compared to microsatellites, and an assignment success
of 97% to the population of origin and of 100% to the region of origin was achieved
when all SNP loci were used. Using outlier SNPs only resulted in minor improve-
ments to assignment success to the population of origin but improved regional assign-
ment. We discuss the implications of these new genetic resources for the conservation
and management of Atlantic salmon in North America.
Keywords: assignment tests, fishery, hierarchical population structure, local adaptation, micro-
satellites, single nucleotide polymorphisms, SNP array
Received 5 July 2014; revision received 14 October 2014; accepted 15 October 2014
Introduction
The increased availability of genomic data in nonmodel
organisms has the potential to revolutionize the use of
genetic tools in conservation and management (Kohn
Correspondence: Jean-S
ebastien Moore, Fax: +1 418 656-7176;
E-mail: jean-sebastien.moore.1@ulaval.ca
1
Co-first authors.
© 2014 John Wiley & Sons Ltd
Molecular Ecology (2014) 23, 5680–5697 doi: 10.1111/mec.12972
et al. 2006; Avise 2009; Primmer 2009; Allendorf et al.
2010). At the most basic level, the greater number of
loci available can increase the precision and accuracy of
population parameters of interest formerly estimated
with a small number of neutral markers (Allendorf et al.
2010). More fundamentally, however, genomic tools can
also provide new types of information not available
with traditional markers. For instance, they can provide
information on the genomic architecture of inbreeding
and/or outbreeding depression (Allendorf et al. 2010)
and highlight regions of differential introgression in the
genome of native taxa under threat from hybridization
(Crispo et al. 2011; Lamaze et al. 2012; Hohenlohe et al.
2013). Furthermore, genomic data can be used to infer
the presence of loci or genomic regions that show
higher differentiation and are putatively under direc-
tional selection (Nosil et al. 2009; Stapley et al. 2010).
Such outlier loci could therefore offer an opportunity to
assess patterns of local adaptation among populations,
thus informing conservation and management.
Establishing conservation units based on objective
and practical criteria has always been a challenge (Mo-
ritz 1994; Waples 1995; Fraser & Bernatchez 2001).
While neutral markers may provide appropriate tools to
delimit reproductively isolated populations (Waples &
Gaggiotti 2006), their utility in determining their adap-
tive potential or ecological distinctiveness has been
questioned (Crandall et al. 2000; Fraser & Bernatchez
2001; McKay & Latta 2002). The inclusion of phenotypic
or ecological data in the definition of conservation units
has been advocated by some (Waples 1991; Crandall
et al. 2000), but remained impractical for most taxa of
conservation concern. The analysis of genomic data
combined with genome scan approaches identifying
putatively adaptive loci could provide a practical solu-
tion to this problem (Bradbury et al. 2013). Recently,
Funk et al. (2012) proposed a framework to integrate
adaptive genetic variation from genomic data into the
definition of conservation units. In short, they propose
a sequential assessment of population structure using in
turn putatively neutral loci and putatively adaptive loci.
In principle, the inclusion of information on putatively
adapted loci can help resolve questions regarding the
spatial scale of local adaptation and help tailor conser-
vation and management actions to specific conditions.
Anadromous Atlantic salmon (Salmo salar) is a species
of biological as well as cultural and economic impor-
tance (Dodson et al. 1998; COSEWIC 2011). Population
declines, however, have also made it a species of con-
servation concern (COSEWIC 2011; Chaput 2012). His-
torically, anthropogenic habitat destruction and
overexploitation have severely reduced the range of the
species in the United States (Spidle et al. 2001) and
reduced abundance in southern Canada (COSEWIC
2011), with more recent declines (since the early 1990’s)
possibly resulting from increased mortality at sea
(Friedland et al. 2005; Chaput 2012). In addition, the
species is popular with the sport fishing industry,
which targets adults returning from their marine migra-
tions to spawn. Offshore fisheries also operate in Green-
land, Labrador and on Saint-Pierre and Miquelon
Islands targeting migrating adults in the marine phase
of their life cycle (Chaput 2012). In Greenland, for
example, a subsistence fishery targets adults migrating
from all regions of North America to the feeding
grounds of West Greenland, making it an important
mixed-stock fishery (Chaput et al. 2005; Gauthier-Ouel-
let et al. 2009). Improving genetic tools for mixed-stock
fishery analysis and individual assignments to natal riv-
ers or to their management unit of origin is therefore of
interest. Furthermore, sampling as many source popula-
tions as possible is crucial to increase the precision of
estimates of stock composition and of individual assign-
ment. The use of a large number of SNPs, and espe-
cially ‘outlier’ SNPs, offers the promise to improve the
power of mixed-stock fishery analysis and assignment
tests (Russello et al. 2011).
An important paradigm in existing conservation and
management strategies for Atlantic salmon is that popu-
lations are structured hierarchically (Dodson et al. 1998;
COSEWIC 2011). At a small scale, the precise homing
behaviour of the species (Stabell 1984; Hendry et al.
2004) leads to genetic differentiation among rivers (King
et al. 2001; Spidle et al. 2001; Dionne et al. 2008; Perrier
et al. 2011; Bradbury et al. 2014a) and in some cases
within rivers (Garant et al. 2000; V
€
ah
€
a et al. 2007; Dion-
ne et al. 2009). At a larger scale, dispersal occurs, but
effective migration is spatially restricted, leading to
genetically differentiated regional groups (Dionne et al.
2008; Tonteri et al. 2009; Bourret et al.
2013a; Bradbury
et al. 2014a). These genetically differentiated regional
groups form the basis of management and conservation
units (COSEWIC 2011), because they represent units
with restricted gene flow that are at least partially
demographically independent (Waples & Gaggiotti
2006). Furthermore, there is increasing evidence that
these genetic groups are differentially adapted to spa-
tially varying environmental conditions (Taylor 1991;
Garcia de Leaniz et al. 2007; Fraser et al. 2011), such as
pathogen diversity (Dionne et al. 2007), water pH (Fra-
ser et al. 2008), temperature and geology (Bourret et al.
2013a). Anadromous Atlantic salmon would thus repre-
sent a good candidate species to evaluate the potential
of adaptive markers to enhance the definition of conser-
vation units.
In this study, data from a geographically extensive
microsatellite database are combined with thousands
of SNP markers to apply a conservation genomics
© 2014 John Wiley & Sons Ltd
CONSERVATION GENOMICS OF ATLANTIC SALMON 5681
framework to anadromous populations of Atlantic sal-
mon in North America. Our first goal was to provide
a North American-wide assessment of population
structure on the basis of both neutral and putatively
adaptive genetic variation to inform the establishment
of conservation units for management. Population
structure has been examined previously for the Qu
ebec
populations (Dionne et al. 2008), the Newfoundland
populations (Bradbury et al. 2014a) and the Nova Sco-
tia populations (McConnell et al. 1997; Vandersteen
Tymchuk et al. 2010). A population genomics assess-
ment conducted on the entire North American range
in a single analysis, however, is lacking (but see King
et al. (2001) for a range-wide assessment of microsatel-
lite DNA variation based on a small number of sam-
pling locations and Verspoor (2005) for a North
American-wide study based on allozymes). By combin-
ing results from previous studies to new data, we now
have the most extensive Atlantic salmon microsatellite
database to date to infer population structure with
high geographical resolution: the database is comprised
of data from 149 sampling locations, with over 9000
individuals genotyped at 15 microsatellite loci. To eval-
uate the usefulness of integrating adaptive variation in
the definition of conservation units, we used SNP data
generated with a medium-density SNP array from a
subset of the populations (Bourret et al. 2013b). Popula-
tion structure was assessed on both putatively adap-
tive and neutral markers identified with two genome
scan approaches. Our second goal was to test the
power of the database for assigning individuals to (i)
regional groups defined using the microsatellite data
and (ii) their sampling location of origin. We did so
with microsatellite markers and SNP markers sepa-
rately, with the explicit aim of evaluating the relative
performance of the two types of markers. For the SNP
data, we also compared the performance of outliers vs.
putative neutral markers in population assignment.
Materials and methods
Microsatellite genotyping and analysis
A total of 9142 anadromous individuals from 149 sam-
pling locations from the entire North American range
of the species were genotyped at 15 microsatellite mark-
ers (Fig. S1 and Table S1, Supporting information). Life
stages sampled include parr, smolts, and in most cases,
returning adults and the samples were collected
between 2000 and 2010 (Table S1, Supporting informa-
tion). The data represented a combination of previously
analysed data sets (see Bradbury et al. (2014a) and Di-
onne et al. (2008) for regional analyses and further
details of genotyping) and new data (Table S1, Support-
ing information). Note that some populations included
have been supplemented by hatchery-reared fish and/
or have been hybridized with farmed escapees (Carr
et al. 2004; Bourret et al. 2011). While this will necessar-
ily influence our results, our goal is to document pres-
ent-day population structure regardless of its causes.
Individuals were genotyped using an ABI 3130xl (or
standardized from ABI 3100 following Gauthier-Ouellet
et al. 2009) by three independent laboratories (Table S2,
Supporting information). Allele scoring was standard-
ized across laboratories in the Laval University labora-
tory using a panel of 10 standard individuals.
Additional individuals (n = 87 total; 46 individuals
from Newfoundland and Labrador; 41 from Nova Sco-
tia and New Brunswick) were genotyped to account for
unrepresented allele values and assess the entire range
of allele sizes. To account for discrepancies between
ABI 3100 and ABI 3130xl, a panel of 64 individuals rep-
resenting all allele sizes of Qu
ebec populations previ-
ously genotyped on the ABI 3100 were analysed on the
ABI 3130xl. Rescreening of samples or re-analysis of
allele sizes was completed when discrepancies existed
between laboratories or machines. Scoring patterns
among laboratories were generally consistent and
allowed standardization using simple rules (see Table
S2, Supporting information).
We tested for departures from Hardy–Weinberg
equilibrium (HWE) using the least-square method
based on the AMOVA F
IS
implemented in GENODIVE
(Meirmans & van Tienderen 2004) with 999 permuta-
tions. We tested for linkage disequilibrium (LD)
between each pair of loci in each population using
Genepop (1000 iterations per batch for 100 batches;
Raymond & Rousset 1995). We adjusted the P-values
of the LD and HWE tests for multiple comparisons
according to the false discovery rate (FDR) method in
the function ‘p.adjust’ in R with an experiment-wide
alpha of 0.05 (R Development Core Team 2006).
Observed (H
O
) and expected heterozygosity (Hs; Nei
1987) were calculated in GENODIVE (Meirmans & van
Tienderen 2004). Allelic richness was calculated for
each sampling location using the rarefaction method
implemented in HP rare (Kalinowski 2005). The Moisie
River sample was removed from this analysis because
it is missing data at three loci. The HP-rare analysis
was repeated twice: once on all the sampling locations
(except Moisie; minimal sample size of 16) and once
with the KC6 and Georges Rivers samples removed
because these two locations had small sample sizes
and many missing values, respectively (minimal sam-
ple size of 36). Pairwise F
ST
values were calculated in
GENODIVE, significance was assessed using 999 per-
mutations, and P-values were adjusted for multiple
comparisons with the FDR method.
© 2014 John Wiley & Sons Ltd
5682 J.- S. MOORE ET AL.
SNP-array genotyping, filters and basic statistics
A total of 1080 anadromous individuals from 50 sam-
pling locations across the North American range (Table
S3, Supporting information) were genotyped at 5568
loci using the SNP array developed by the Centre for
Integrative Genetics (CIGENE, Norway) following the
manufacturer’s instructions (Illumina, San Diego, CA,
USA). The data from the Qu
ebec populations were pre-
viously published in (Bourret et al. 2013a), but data
from the other regions are previously unpublished.
Detailed methods for SNP discovery and quality control
can be found in Bourret et al. (2013a). Ascertainment
bias for the North American populations is minimal
(Bourret et al. 2013b) and should not bias our results
because most markers on the chip were discovered
from populations not considered here (European). The
quality of individual samples was assured by only
using individuals genotyped at a >0.95 call rate (CR:
proportion of genotyped SNPs). Markers with minor
allele frequencies less than 1% (MAF < 0.01) across all
populations and markers missing in more than 5% of
individuals were excluded from our analyses. Observed
(H
O
) and expected heterozygosity (H
S
; Nei (1987)) were
calculated in GENODIVE (Meirmans & van Tienderen
2004). As with the microsatellite data, we tested for
departures from HWE GENODIVE with 999 permuta-
tions and adjusted the P-values for FDR in R. Pairwise
F
ST
values were calculated in GENODIVE, and signifi-
cance was assessed using 999 permutations.
Outlier markers detection
We used several alternative genome scan methods to
detect loci with greater than expected levels of diver-
gence among regional groups. First, we used BAYE-
SCAN, a Bayesian approach that allows the estimation
of the posterior probability of a given locus being under
the effect of selection (Foll et al. 2008). It is based on the
multinomial-Dirichlet model, and assumes that allele
frequencies among demes are correlated through a com-
mon migrant gene pool, therefore allowing complex
ecological scenarios to be modelled satisfactorily (Foll
et al. 2008). Because previous studies suggested that
local adaptation is probably important at the regional
scale (Dionne et al. 2008; Bourret et al. 2013a), we ran
BAYESCAN on the entire data set by defining groups
of individuals on the basis of the microsatellite-defined
regional groups (not on populations; Table S1, Support-
ing information and results) using all the defaults (tests
runs with longer chain parameters gave identical
results). Loci putatively under selection were defined as
those with alpha-values significantly >0 (i.e. with Q-val-
ues smaller than 0.05) while loci putatively under
balancing selection had alpha-values significantly smal-
ler than 0. All other loci were considered neutral. Sec-
ond, we used hierarchical Fdist (Excoffier et al. 2009)
implemented in Arlequin 3.5 (Excoffier & Lischer 2010).
This method uses a hierarchical island model allowing
for lower migration rates among groups compared to
among populations within groups. This has been shown
in some instances to reduce the prevalence of false posi-
tives compared to methods using the finite island
model (Excoffier et al. 2009; Mita et al. 2013). We used
the regional groups defined using the microsatellites
(Table S1, Supporting information and results) as the
higher level of population structure. Third, as an alter-
native way to identify loci that maximize assignment to
the sampling location of origin, we conducted a ‘nested
Fdist’ analysis, where we used a nonhierarchical Fdist
on each microsatellite-defined regional groups sepa-
rately, and combined all the unique markers identified
by each of the individual analyses. The goal with the
nested analysis was to select markers that would maxi-
mize assignment power to the population and is thus
only discussed in this context. In both hierarchical and
nested Fdist analyses, loci with significantly (at the 0.01
significance level) higher F
CT
or F
ST
values were classi-
fied as outliers potentially under directional selection
among regional groups, loci with significantly lower
F
CT
or F
ST
values were loci putatively under balancing
selection, and all other loci were classified putatively
neutral loci. The results of both the Bayescan and hier-
archical Fdist analyses were plotted against the esti-
mated genomic position of a subset of SNP loci for
which map data were available (Brenna-Hansen et al.
2012) to assess whether specific genomic regions
or linkage groups showed elevated amounts of
divergence.
Population structure
Initial tests using the widely used Bayesian clustering
program STRUCTURE (Pritchard et al. 2000) deter-
mined that this approach was too computationally
demanding even for the microsatellite data set. Instead,
we used BAPS (Corander et al. 2003, 2004, 2008), a
model-based Bayesian clustering approach that, like
STRUCTURE, infers genetic groupings that maximize
Hardy–Weinberg and linkage equilibrium (Corander
et al. 2003). Unlike STRUCTURE, however, BAPS infers
the optimal number of clusters directly (instead of rely-
ing on ad hoc measures; Corander et al. 2004) and is
computationally more efficient, therefore facilitating the
analysis of large data sets. Because in this analysis we
were interested in regional clustering of sampling loca-
tions, we used the ‘clustering of groups of individuals’
© 2014 John Wiley & Sons Ltd
CONSERVATION GENOMICS OF ATLANTIC SALMON 5683
option in BAPS. This prior is justified by the homing
behaviour of Atlantic salmon (Stabell 1984; Hendry
et al. 2004). For the microsatellite data set, 20 replicate
runs of BAPS were performed with a maximum num-
ber of K of 149 (the number of samples in the data set).
Runs with smaller values for maximum K (10, 20, 30, 40
and 50) were also performed to explore the effects
of this prior on the final results, which were always
identical.
After running BAPS on the entire data set, we exam-
ined the possibility of further hierarchical population
structure, by running both BAPS and STRUCTURE on
the regional groups defined by BAPS with the full data
set. For these analyses on each regional groups, STRUC-
TURE was run with the admixture model with corre-
lated allele frequencies, 50 000 burn-in and 100 000
MCMC repetitions, K 1–10 with 10 iterations of each K.
Alpha-values were examined to ensure that conver-
gence had been attained after the burn-in (Pritchard
et al. 2000). For the BAPS analysis, we ran 20 indepen-
dent runs with max K = 50 on each of the regional
groupings separately.
For the SNP data set, BAPS was run on all SNPs, on
the neutral SNPs identified by Bayescan, on the outlier
SNPs identified by Bayescan and on the outlier SNPs
identified by Fdist. In all cases, 20 independent runs
were performed with a maximum number of K of 50
(again, the number of sampling locations).
Finally, to infer the relationship of regional groups
among each other, we constructed an unrooted neigh-
bour-joining tree of Cavalli-Sforza chord distances for
the microsatellite data using the software package Phy-
lip (Felsenstein 1993). The Moisie was excluded from
this analysis because it is missing data at three loci. We
used 1000 bootstrap replicates to assess confidence in
the nodes.
Individual assignment to regional groups and river of
origin
We used two approaches to test and compare the rela-
tive power of the microsatellite and SNP data sets to
assign individuals to their region of origin and their
sampling location of origin. For the definition of regio-
nal groups, we used two hierarchical levels of popula-
tion structure defined by the microsatellites (see
Results). First, we used the leave-one-out approach
implemented in ONCOR (Kalinowski et al. 2007). This
approach sequentially removes each individual from the
baseline and estimates its origin by determining the
sampling location that has the highest probability of
producing its genotype. Assignment success is then esti-
mated as the proportion of individuals that were cor-
rectly assigned to their sampling location of origin. This
approach does not allow for missing data. We therefore
filled the missing data by drawing alleles randomly
from the entire data set in GENODIVE (Meirmans &
van Tienderen 2004). This will reduce overall population
differentiation and will in most cases provide conserva-
tive estimates of power of assignment. Second, we used
the assignment test option in GENODIVE, which imple-
ments the approach of Paetkau et al. (1995, 2004). The
assignment tests were performed on two different mi-
crosatellite data sets and on four different SNP data sets.
For the microsatellites, analyses were first conducted on
the entire database. However, to minimize biases due to
sampling effects in the comparisons between SNPs and
microsatellites, we combined the different sampling
locations that represented several tributaries of the same
river (see Supporting information for details). The sec-
ond microsatellite data set tested only included the sam-
pling locations for which we had SNP data. This was
performed to further ensure comparability between the
microsatellite and SNP assignment tests in case the
inclusion of many genetically similar populations in the
full microsatellite data set biased assignment success.
For the SNPs, our goal was to compare the power of the
entire data set to other data sets, which include only
outlier loci and which have been found elsewhere to be
sufficient for powerful assignment (Russello et al. 2011).
We therefore performed assignment tests on a data set
including all loci (n = 3192), using only the outliers from
the Bayescan analysis (n = 106), the outliers from the
hierarchical Fdist analysis (n = 61) and from the nested
hierarchical Fdist (n = 293). We further evaluated the
effect of the number of SNPs included in the analysis on
the power of assignment tests. We tested data sets with
50, 100, 250, 500 and 1000 randomly selected SNP mark-
ers and performed assignment tests in ONCOR for three
replicate data sets for each number of random SNPs at
both the population and regional level. In all cases,
missing values were filled in GENODIVE according to
the overall allele frequencies. Note that using the same
samples to identify outliers and to assess their power of
assignment can results in upwardly biased estimates of
power (Anderson 2010). Our results, however, show that
adding those outliers did not improve our assignment
success compared to using all other markers (see
Results), and the potential upward bias in power only
reinforces this conclusion.
Results
Microsatellites
The microsatellite loci were highly polymorphic, with
the number of alleles per locus ranging from 14 (D486)
to 91 (SsaD144) (average 40; total 593). Missing values
© 2014 John Wiley & Sons Ltd
5684 J.- S. MOORE ET AL.
