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).