Sample collection and processing
All starlings sampled were collected during the non-breeding season,
when large flocks aggregate at high-quality food sources. These samples
therefore do not represent discrete breeding populations but rather a
mixture of birds from the surrounding region, potentially including
migrants from more northerly areas. Starlings (N = 166) were
collected in January-March 2016 and 2017 from 26 dairies and feedlots by
the U.S. Department of Agriculture’s Wildlife Services personnel in
Arizona, California, Colorado, Idaho, Illinois, Iowa, Kansas, Missouri,
Nebraska, Nevada, New Hampshire, New Mexico, New York, North Carolina,
Texas, Washington, and Wisconsin. USDA Wildlife Services personnel
collected birds from 2-5 sites in each state, where each collection site
was >5 km apart. USDA personnel then euthanized whole adult
males by lethal gunshot, avicide, or live traps, stored at these birds
at 0°C until tissue sampling, and recorded the latitude and longitude of
each collection site using a handheld GPS. The collection and use of
starlings for this and related studies were approved by the U.S.
Department of Agriculture, National Wildlife Research Center’s
Institutional Animal Care and Use Committee (QA-2572, S.J. Werner -
Study Director; Werner, Fischer, & Hobson, 2020; Table S1).
Breast muscle tissue was sampled using biopsy punches (Integra Miltex,
York, PA) and frozen in 95% ethanol. Samples were shipped on dry ice,
and DNA was extracted using a Qiagen DNeasy kit following the
manufacturer’s protocol (Qiagen, New York, NY). DNA concentration of
each sample was quantified using a Qubit 2.0 fluorometer (Thermo Fisher
Scientific, New York, NY). Following the protocol of Peterson et
al. (2012), we generated a reduced-representation genomic dataset of
double-digested, restriction-site associated DNA markers as described in
Thrasher et al. (2017) using the restriction enzymes SbfI and
MspI and adaptors P1 and P2. We sequenced 100bp, single-end reads on an
Illumina HiSeq 2500. We trimmed and filtered for quality using the
FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit). We then used
the process_radtags commands in stacks v 1.19 (Catchenet al., 2013) to demultiplex the remaining sequences. In
subsequent filtering steps, we retained reads only if the following
conditions were met: reads passed the Illumina chastity filter,
contained an intact Sbf I RAD site, contained one of the unique
barcodes, and did not contain Illumina indexing adaptors.
We assembled sequences to an S. vulgaris reference genome
(Hofmeister, Rollins et al., in prep ) using the ref-map
option in stacks. Individual reads were mapped to the reference
genome using bowtie2 version 2.2.8 (Langmead & Salzberg, 2012)
using the “very sensitive local” set of alignment pre-sets. The
bioinformatics pipeline used for the reference-based assembly has the
advantage of using fewer similarity thresholds to build loci. We
required that a SNP be present in a minimum of 80% of the individuals
(-r 0.8) with a minimum stack depth of 10 reads at a locus within an
individual (-m 10) for it to be called. We removed two individuals, one
with >50% missing data and one with >50%
relatedness (measured using the unadjusted AJK statistic and calculated
within vcftools ), leaving 158 individuals remaining in the study.
We also used vcftools –hwe option to remove any SNPs out of
Hardy-Weinberg equilibrium (e.g., an exact test that compared expected
and observed heterozygosity gave a P-value less than 0.001). About 6%
of sequenced variants were out of HWE, and this filter left 14,134
variants remaining. For all genotype-environment analyses we used all
SNPs in a given stack, but for STRUCTURE and other analyses sensitive to
linkage disequilibrium we used only the first SNP in each stack; for
unphased genomic data like the RAD loci analyzed here, this strategy of
exporting one SNP per stack is often used as an indirect method of
controlling for linkage.
Patterns of genetic diversity and differentiation
We estimated per-locus measures of genetic diversity and genome-wide
differentiation using the populations option in Stacks (Catchenet al ., 2013). We used vcftools to calculate
FST among population pairs and heterozygosity and
nucleotide diversity (π) within populations (Danecek et al.,2011). We determined expected heterozygosity at the population level
using the Hs() function in the R package adegenet (Jombart
,2008), and tested for differences in observed and expected
heterozygosity of each locus using a paired t-test in base R (R Core
Team, 2019). We investigated genetic structure within North American
starlings using an analysis of molecular variance (AMOVA) in the R
package poppr (Kamvar et al., 2014). We tested whether
most genetic variation was observed among individuals or among sampling
sites (“populations”). To determine significance, we compared observed
variation at each hierarchical level to the randomly permuted distance
matrices for that particular level using 1000 simulations in the
function randtest() in the R package adegenet (Jombart,
2008), hypothesizing that the observed variance is greater than expected
within individuals and less between individuals and between populations.
We tested for isolation by distance (IBD) using a simple Mantel test inadegenet (Jombart, 2008): for these data, the assumption of
stationarity likely holds, given that North American starlings appear to
be in mutation-migration-drift equilibrium (Guillot & Rousset, 2013).
Geographic distances among sampling locations follow a bimodal
distribution where locations are either very near or far from each
other, and thus we caution that the Mantel test may not be an
appropriate test of isolation by distance (Supplementary Information).
Nevertheless, we report the results of a Mantel test with a Spearman
correlation and 10000 permutations. Finally, we test for isolation by
environment (IBE) using the R package vegan to compare
environmental, geographic, and genetic distances (Oksanen et al.,2019). Controlling for geographic distance, we ran partial Mantel
correlograms using Euclidean environmental distance matrices built from
all bioclimatic variables retained in the Selection analysesdescribed in Section 4 below.
Population structure
We first used non-parametric approaches to determine whether starlings
clustered by sampling location, using principal components analysis inSNPRelate (Zheng et al. , 2012) and discriminant analysis
of principal components (DAPC) in adegenet (Jombart, 2008). We
then tested for population structure using STRUCTURE (Pritchard et
al., 2000) by simulating 10 runs of each K . Although we sampled
starlings from 17 locations, we hypothesized that North American
starlings would cluster into at most eight populations (K =1-8)
based on the non-parametric tests described above. To select the
best-supported K, we used the Evanno method implemented in STRUCTURE
HARVESTER v0.6.94 (Earl & vonHoldt, 2011). We averaged results across
the 10 runs using the greedy algorithm in the program CLUMPP v1.1.2
(Jakobsson & Rosenberg, 2007), and visualized results using DISTRUCT
v1.1 (Rosenberg, 2003). Given that evidence of population structure
depends on the filtering thresholds selected, we ran this model-based
approach using a very strict minimum minor allele frequencies (MAF=0.3)
and a more relaxed minimum frequency (MAF=0.1) (Linck & Battey, 2017).
STRUCTURE results did not differ substantially among MAF thresholds, and
we used a filtered dataset with a minimum minor allele frequency of 0.1
in the demographic analyses below.
Demographic modeling
We used a traditional model-based approach (fastsimcoal2) to explicitly
test for a bottleneck in North American starlings (Excoffier et
al., 2013). We first estimated the folded SFS using a Python script
from Simon Martin (available at
https://github.com/simonhmartin/genomics_general). We modified the
simple 1PopBot20Mb.tpl and .est files provided with the software as
follows: the size of the bottleneck (NBOT: 10-1000), start of bottleneck
(TBOT: 10-300), ancestral size (ANC: 1000-100000), current size of the
population (NCUR: 100-100000), and end of bottleneck (TENDBOT: TBOT +
500). We performed 100 runs of this model and assessed model fit using
delta-likelihood for each run of the model; since all models use the
same numbers of parameters, the Akaike information criteria do not
change between models. We also ran a demographic model that takes
linkage into account, the Stairway plot method, but given that this
method does not explicitly test for a genetic bottleneck, we describe
the method and results in the Supplementary Information.
Selection analyses
We first identified environmental variation at each sampling location
using the R package raster (Hijmans & van Etten, 2012),
extracting all 19 bioclimatic variables for each set of sampling
coordinates from WorldClim 2.0 at a resolution of 5 min on June 16, 2018
(Fick & Hijmans, 2017). We recognize that these environmental
conditions represent those experienced during the non-breeding season,
and therefore they do not represent the full range of conditions
experienced by the starlings sampled; see the Supplementary Information
for a discussion of our assumptions about this environmental sampling.
We used redundancy analysis (RDA) to examine how loci may covary across
multiple environmental gradients (Forester et al., 2018). RDA is
especially powerful when testing for weak selection, and detects true
positives in large data sets more reliably than other multivariate
methods like Random Forest (Forester et al., 2018). Because RDA
requires no missing data, we first imputed genotypes where missing sites
were assigned the genotype of highest probability across all
individuals—a conservative but quick imputation method. Before
imputation and after filtering, about 31% of genotypes were missing
from the dataset of 15,038 SNPs (not filtered for MAF), and we required
that all imputed SNPs were present in at least 80% of individuals. We
then used the R package vegan to run the RDA model (Oksanenet al. , 2019); for a full description of this method, see
Forester et al., 2018. Briefly, RDA uses constrained ordination
to model a set of predictor variables, and unconstrained ordination axes
to model the response (genetic variation). RDA infers selection on a
particular locus when it loads heavily onto one or more unconstrained
predictor axes, and we retained five variables with relatively low
variance inflation factors: BIO1 or Annual Mean Temperature (VIF=3.54),
BIO7 or Temperature Annual Range (4.55), BIO12 or Annual Precipitation
(8.69), BIO16 or Precipitation of Wettest Quarter (7.91), and elevation
(2.19). We tested for significance using the anova.cca function
within the vegan package, and also permuted predictor values
across individuals to further check significance of the model. We
identify candidate loci as those that are 3 or more standard deviations
outside the mean loading. The R scripts for all RDA analyses and figures
were written by Brenna Forester (available at
https://popgen.nescent.org/2018-03-27_RDA_GEA.html).