Habitat-linked genetic variation supports microgeographic adaptive divergence in an island-endemic bird species.

We investigated the potential mechanisms driving habitat-linked genetic divergence within a bird species endemic to a single 250-km2 island. The island scrub-jay (Aphelocoma insularis) exhibits microgeographic divergence in bill morphology across pine-oak ecotones on Santa Cruz Island, California (USA), similar to adaptive differences described in mainland congeners over much larger geographic scales. To test whether individuals exhibit genetic differentiation related to habitat type and divergence in bill length, we genotyped over 3000 single nucleotide polymorphisms in 123 adult island scrub-jay males from across Santa Cruz Island using restriction site-associated DNA sequencing. Neutral landscape genomic analyses revealed that genome-wide genetic differentiation was primarily related to geographic distance and differences in habitat composition. We also found 168 putatively adaptive loci associated with habitat type using multivariate redundancy analysis while controlling for spatial effects. Finally, two genome-wide association analyses revealed a polygenic basis to variation in bill length with multiple loci detected in or near genes known to affect bill morphology in other birds. Our findings support the hypothesis that divergent selection at microgeographic scales can cause adaptive divergence in the presence of ongoing gene flow.

INTRODUCTION

Understanding the processes that generate and maintain adaptive phenotypic variation is a fundamental goal of evolutionary biology. Much research has focused on the interaction between gene flow and selection in shaping patterns of local adaptation across selective landscapes (Haldane, 1930; Räsänen & Hendry, 2008; Tigano & Friesen, 2016). The homogenizing effects of gene flow on local adaptation (Akerman & Bürger, 2014; Slatkin, 1973, 1987) have led to the dominant paradigm that isolation is generally a prerequisite for adaptive divergence among populations (Garant et al., 2007; Kawecki & Ebert, 2004; Lenormand, 2002; Nosil, 2008; Yeaman & Whitlock, 2011). Yet there is mounting evidence that “microgeographic” adaptation in the absence of geographic barriers may be more common than predicted by traditional models (reviewed by Richardson et al., 2014).

Microgeographic adaptation is a specific case of local adaptation that occurs within the “dispersal neighborhood” of a species, defined as two standard deviations of the dispersal distribution of a population (Richardson et al., 2014; Wright, 1943, 1946). Early examples of microgeographic adaptation involved strong divergent selection over very small spatial scales in plants (Antonovics, 2006; Turner et al., 2010). Recent evidence indicates that adaptive divergence driven by strong divergent selection is also possible in highly mobile animals with significant levels of dispersal and gene flow (e.g., Hohenlohe et al., 2010; Mikles et al., 2020; Nacci et al., 2016; Torres‐Dowdall et al., 2012; Urban et al., 2017). While many studies clearly demonstrate phenotypic and genetic variation consistent with hypotheses of microgeographic adaptive evolution (e.g., Charmantier et al., 2016; Maciejewski et al., 2020; Pequeno et al., 2021), determining the environmental factors and evolutionary and genetic mechanisms driving these patterns remains a difficult challenge (Barrett & Hoekstra, 2011; Hoban et al., 2016). Studies that are successful at showing both genomic evidence of divergent selection and a genetic basis to diverging phenotypes at fine spatial scales are generally restricted to traits controlled by few genes of large effect (e.g., Laurent et al., 2016; Linnen et al., 2013; Nosil et al., 2018; Pfeifer et al., 2018). However, many traits under selection are polygenic, making multiple loci of small effect difficult to detect in common selection tests such as F ST outlier methods (Hoban et al., 2016; Lundregan et al., 2018; Tiffin & Ross‐Ibarra, 2014; Wellenreuther & Hansson, 2016; Yeaman, 2015).

Spatial population genetic structure is another complicating factor in studies of local adaptation (Battey et al., 2020). Theoretical and empirical research has demonstrated that gene flow can erode genetic signatures of local adaptation by replacing locally adapted alleles with maladaptive immigrant alleles (Hendry et al., 2002; Kawecki & Ebert, 2004; Lenormand, 2002; Postma & van Noordwijk, 2005; Räsänen & Hendry, 2008). Recent work, however, has illustrated how nonrandom dispersal with respect to genotype has the potential to mitigate maladaptive gene flow (Edelaar & Bolnick, 2012; Edelaar et al., 2008; Lowe & McPeek, 2014; Nicolaus & Edelaar, 2018). Nonrandom dispersal may involve multiple mechanisms, including individuals avoiding habitats to which they are not adapted (Edelaar et al., 2019; Wang & Bradburd, 2014), reduced fitness in nonoptimal habitat (Karpestam et al., 2012; Richardson & Urban, 2013), assortative mating (Servedio, 2016) or a combination of mechanisms (e.g., “magic trait,” where a trait undergoing divergent selection also contributes to assortative mating; Servedio et al., 2011). We can therefore predict that nonrandom dispersal and subsequent directional gene flow would facilitate adaptation due to increased genome‐wide differentiation (Edelaar & Bolnick, 2012). Such segregation of adaptive genetic variation has even been demonstrated across fine spatial scales in mobile taxa (e.g., Bolnick et al., 2009; Camacho et al., 2013, 2020), and could be a potential mechanism facilitating microgeographic adaptation in some systems.

While an increasing number of studies provide evidence for microgeographic adaptation (Richardson et al., 2014), conclusively demonstrating fine‐scale adaptation requires a highly integrative analytical framework that provides multiple lines of evidence. First, adaptive trait variation must be documented across environmental gradients and within the dispersal neighbourhood of the species in question. This background knowledge of the study species also provides a priori hypotheses about which environmental variables are most important for local adaptation. Second, there must be evidence for divergent selection across the environmental gradient that could drive the observed trait divergence (Ahrens et al., 2018; Blanquart et al., 2013; Richardson & Urban, 2013; de Villemereuil et al., 2014). Last, some of the observed phenotypic variation must be genetically based (Hendry et al., 2001; Reid et al., 2016; Turner et al., 2010). Recent studies have demonstrated that complex, multilocus patterns of selection can be detected in heterogeneous landscapes in the presence of gene flow by using genome scan approaches (Forester et al., 2016, 2018; Lotterhos & Whitlock, 2015; Yeaman, 2015). We use such a framework in this study.

The island scrub‐jay (Aphelocoma insularis), a medium‐sized bird restricted to Santa Cruz Island, California, USA, may represent a case of microgeographic adaptation. This species has been evolving in isolation for ~1 million years (Delaney & Wayne, 2005; McCormack et al., 2011) with no evidence of gene flow from outside the island to confound genetic patterns of differentiation, as seen in other island systems (Postma & van Noordwijk, 2005; Postma et al., 2009). The 250‐km2 island is variable in topography, ecology and climate. Coniferous forest was predominant on the island during the Pleistocene (Anderson et al., 2009). Current vegetation cover on the island is a mosaic of coastal sage scrub, oak woodlands, and chaparral habitat dominated by island scrub oak (Quercus pacifica; Junak, 1995). Only three relict stands of bishop pine (Pinus muricata) remain in the western, central and eastern regions of the island (Fischer et al., 2009; Walter & Taha, 1999; Figure 1).

FIGURE 1

Study area of Santa Cruz Island, California, USA, with sampling sites (black dots; n = 152) and habitat type (green, pine; orange, oak). Island scrub‐jays were sampled from each of the pine–oak ecotones (western oak = 39, central oak = 60, eastern oak = 12, western pines = 25, central pines = 7, eastern pines = 9). Insets show regional breakdown of Santa Cruz Island and location of the island (in black) relative to mainland California

Despite its very restricted range, the island scrub‐jay exhibits phenotypic divergence in bill size and shape between different habitat types. Jays that occur in oak‐dominated habitat have shorter, deeper bills compared to jays located in pine habitat (Langin et al., 2015). This striking pattern mirrors well‐described adaptations seen in mainland populations of California (A. californica) and Woodhouse's scrub‐jays (A. woodhouseii; Peterson, 1993), which are thought to be adaptive, as the short, stout bills of jays living in oak woodlands are more efficient for hammering open acorns, while long, shallow bills of jays living in pine–juniper forests are more efficient at extracting seed from pine cones (Bardwell et al., 2001). While this divergence is observed in mainland jays separated by hundreds of kilometres, the same pattern occurs across a few kilometers in island scrub‐jays (Langin et al., 2015).

Here we investigate if this pattern of fine‐scale phenotypic divergence in island scrub‐jay bill morphology represents a case of microgeographic adaptation, which would require demonstrating divergent selection across the pine–oak ecotone, and whether the variation in bill morphology has a genetic basis, as seen in other avian systems (Abzhanov et al., 2006; Bosse et al., 2017; Lundregan et al., 2018). Prior studies using pedigree‐based analyses suggested bill length and depth to be heritable (Langin et al., 2015). This previous work used a limited number of microsatellites to elucidate population structure on the island, but was not designed to test for divergent selection or identify loci underlying bill morphology (Langin et al., 2015). Genomic approaches that allow genotyping thousands of markers, such as single nucleotide polymorphisms (SNPs), provide greater statistical power to detect fine‐scale genetic divergence (Allendorf, 2017; Hohenlohe et al., 2021; Morin et al., 2004). Genomic data also allow detection of divergent selection between habitats, to test the role of selection in driving phenotypic divergence (Hoban et al., 2016), in addition to testing whether trait variation is genetically based using genome‐wide association analyses (GWAs; Santure & Garant, 2018).

We genotyped 123 adult male island scrub‐jays using restriction‐site associated DNA sequencing (RADseq) to address three questions: (i) Is gene flow restricted between habitat types? (ii) Does divergent selection act across habitat types? (iii) Does divergence in bill morphology have a genetic basis? We first applied a neutral landscape genomics approach to test alternative hypotheses regarding the effects of variation in habitat type (pine vs. oak) and other landscape features on genome‐wide divergence. We predicted observed genetic divergence would be consistent with limited dispersal or nonrandom gene flow between pine and oak habitats, while also considering the effects of geographic distance and habitat quality. Second, we used a multivariate genotype–environment association (GEA) analysis to test if divergent selection across the pine–oak ecotone drives microgeographic divergence in bill morphology. Finally, we used two different GWA methods to test the hypothesis that variation in bill morphology is genetically based, and that loci correlated with bill morphology are linked to candidate genes that are known to cause variation in bill morphology in other bird species.

METHODS

Island scrub‐jay sampling

Island scrub‐jay blood samples were collected from 152 adult males during September–December, 2009–2011 (see Langin et al., 2015, for detailed field methods; Figure 1). We focused on males to eliminate variation in bill length driven by sex. Each captured jay was measured by the same person (K. Langin) using digital calipers to record: bill length (to ±0.01 mm), measured from the anterior end of the nares to the tip of the bill; bill depth, measured at the anterior end of the nares; and tarsus length. Wing chord and tail lengths were also measured with a ruler (to ±0.5 mm). All work with living birds was approved by the Institutional Animal Care and Use Committees at Colorado State University (IACUC: #887) and the Smithsonian Institution.

RADseq and genotyping

We used RADseq to genotype thousands of anonymous loci throughout the island scrub‐jay genome. We extracted genomic DNA from blood using the Qiagen DNeasy Blood and Tissue Kit (QIAGEN). The samples were individually barcoded and pooled following the protocol from Etter et al. (2011). We generated RADseq libraries using the restriction enzyme SbfI, which targets an 8‐bp cut site (5′ CCTGCAGG 3′). We ligated a unique 6‐bp DNA barcode to each individual's cut DNA before sample multiplexing was performed in equimolar proportions by groups of 14 individuals per pool. We sheared 80–100 µl of each pooled sample to an average fragment size of 500 bp using a Covaris ultrasonicator (Covaris Inc.). Desired fragment sizes between 300 and 600 bp were separated by manual gel excision from the sheared DNA. We then sequenced 100‐bp single‐end reads on two Illumina HiSeq 2500 lanes at the Genomics Core Facility at the University of Oregon (gc3f.uoregon.edu).

Sequenced reads were quality‐filtered and demultiplexed, and individual barcodes were removed using the “process_radtags” program in stacks version 2.3b (Catchen et al., 2013; Rochette & Catchen, 2017; Rochette et al., 2019). We aligned sequences to the Florida scrub‐jay (Aphelocoma coerulescens) draft reference genome (Feng et al., 2020) using bwa‐mem (version 0.7.17; Li & Durbin, 2009). We used the default settings of ref_map.pl in stacks to identify SNPs and exported the raw SNP matrix using the “populations” module of stacks (Catchen et al., 2013). We then visualized and evaluated this matrix for missingness and genotype miscall rates using the R packages genoscapertools (Anderson, 2020a; https://github.com/eriqande/genoscapeRtools) and whoa (Anderson, 2020b; https://github.com/eriqande/whoa). After this evaluation, individuals with >24% missing data were removed by rerunning populations with a whitelist.

Additional filtering was performed using the radiator package (version 1.1.5; Gosselin et al., 2020) in r (R Core Team, 2013). We filtered by global minor allele count (MAC = 3) to reduce sequencing errors, assembly artefacts and rare alleles (Linck & Battey, 2019). We then retained loci with coverage between 6× and 100× across individuals, removing loci with coverage too low for accurate genotype calling and loci with high coverage that probably reflect repetitive regions (e.g., paralogues). We retained markers that were genotyped in at least 70% of individuals and between positions 1 and 88 (removing the last seven positions of our reads, which showed signatures of low‐quality genotyping). We then selected the SNP with the highest minor allele count per contig to minimize effects of short‐distance linkage disequilibrium. SNPs with significant deviation from Hardy–Weinberg proportions (p < .0001) across all samples were removed.

As the Florida scrub‐jay draft genome is a scaffold‐level assembly and unannotated (scaffold N50 = 7.7 Mb), we mapped scaffolds to individual chromosomes of the zebra finch genome assembly (GenBank accession GCA_008822105.2) using the default parameters of satsuma synteny version 2.1.0 (Grabherr et al., 2010). We then used custom R scripts modified from Van Doren et al. (2017) to reorder our island scrub‐jay VCF file relative to the zebra finch genome and remove SNPs where chromosomal positions could not be determined. We imputed missing genotype values for all SNPs (7.01% missing) using beagle (version 5.1; Browning & Browning, 2016) with 25 iterations and an N e setting of 370 (see Section 3). This imputed data set was used for analyses requiring complete data frames: principal components analysis (PCA; Section 2.3), GEA tests (Section 2.4) and genome‐wide associations (GWA; Section 2.5).

Is gene flow restricted between habitat types?

After filtering, we identified and removed loci showing signatures of selection using pcadapt (Luu et al., 2016) with a false discovery rate (FDR) of 10% to avoid confounding neutral demographic patterns with patterns generated by loci under selection when conducting downstream neutral population genetics analyses. We calculated population genomic statistics using 3408 neutral markers across all sampled individuals (Figure 1) analysed as a single population. Observed and expected heterozygosity, nucleotide diversity and the inbreeding coefficient were calculated using the populations program in stacks. Effective population size (N e) was estimated using the linkage disequilibrium method of Waples et al. (2016), calculated in the R package stratag (Archer et al., 2016). We calculated two individual‐based genetic distances: genetic distance based on relatedness (Smouse & Peakall, 1999) using the R package popgenreport (Adamack & Gruber, 2014), and the proportion of shared alleles (Bowcock et al., 1994) using adegenet (Jombart, 2008; Jombart & Ahmed, 2011).

We tested for population genetic structure across the island and between pine and oak habitats using a combination of methods. First, we tested for isolation‐by‐distance using both individual measures of genetic distance outlined above, and pairwise geographic distances calculated using the geodist package to measure “geodist” (Padgham & Sumner, 2020). We tested for isolation‐by‐distance using a Mantel test with 1000 permutations and the “pearson” method in the vegan r package (Oksanen et al., 2013). Second, we used PCA implemented in vegan to visualize how genetic variation is distributed across a reduced number of orthogonal axes without an underlying assumption of genetic groups or spatial structure.

We tested for discrete population structure using snmf in the lear package (Frichot & François, 2015) and structure (Pritchard et al., 2000). Because there are three pine stands, we conducted runs for K = 1–6 for both analyses to test the hypothesis that each pine and adjacent oak stand represent distinct populations. We ran with 10 replicates per K value and retained the run with the smallest cross‐entropy to choose the optimal number of clusters (Frichot et al., 2014). The best supported value of K was the estimate for which the cross‐entropy curve exhibited a plateau or a clear minimum value. We ran structure with a burn‐in length of 50,000 and a run length of 50,000 using the admixture model. We assumed K = 1 if its log likelihood value was highest, as the second‐order rate of change in log likelihood cannot be calculated for K = 1 (Evanno et al., 2005).

Finally, because Santa Cruz Island is characterized by extremely variable and rugged terrain, we predicted that landscape features could influence dispersal and gene flow (Langin, Sillett, Funk, et al.,  2017; Sillett et al., 2012). To understand how different landscape features shape the genetic variation of island scrub‐jays, we applied a linear mixed effects modelling approach to test the effects of landscape features on neutral population structure using maximum likelihood of population effects (MLPE; Clarke et al., 2002; Row et al., 2017; van Strien et al., 2012). MLPE uses individual‐based genetic distances as the response and environmental resistances as fixed effects, with an additional random effect matrix of pairwise individual comparisons to control for their nonindependence (Balkenhol et al., 2016; Trumbo et al., 2019). Although MLPE has received criticism regarding its power to accurately capture how landscape features impact gene flow (Peterman & Pope, 2021), our goal here was to test the validity of a priori hypotheses regarding scrub‐jay dispersal across habitat types (Langin et al., 2015), rather than finding the best supported resistance model. Additionally, MLPE is an effective method to infer associations between landscape features and gene flow (Kozakiewicz et al., 2019; Shirk et al., 2017; Trumbo et al., 2019). We used the proportion of shared alleles and relatedness as our genetic distance response variables. The environmental resistances used as the explanatory variables were those that we hypothesized a priori are related to genetic distance in this system: topographically corrected geographic distance (positive effect), low vegetation density (positive effect) and elevation (positive effect; Figure S1). Landscape data were converted to resistance surfaces using the Reclassify and Raster Calculator tools in arcgis version 10.6 (ESRI). We calculated vegetation density from chlorophyll reflectance in visual and near‐infrared spectra (i.e., enhanced vegetation index) using moderate‐resolution (250 m) imaging spectroradiometer data collected in 2015 (modis.gsfc.nasa.gov). Environmental resistances among individuals were calculated using circuitscape (McRae, 2006) for each landscape resistance surface. To test for resistance between pine and oak habitats, we calculated the percent pine and percent oak within a 300‐m radius of each scrub‐jay sampling location (the diameter of the largest island scrub‐jay territory; Caldwell et al., 2013) using a reclassified 2005 vegetation map of Santa Cruz Island (Langin et al., 2015; The Nature Conservancy, 2007). We then calculated the absolute differences in pine and oak composition by subtracting the percent pine and percent oak between all individuals to help identify transitions between these habitat types (ecotone; positive effect). Because collinearity can cause instability in parameter estimation in regression models (Dormann et al., 2012; Row et al., 2017), we tested each landscape variable for multicollinearity, both prior to and after calculating environmental resistances in circuitscape, to ensure Pearson's r correlations < .7 and variance inflation factor (VIF) scores <5 in our final models.

We standardized landscape resistances to units of standard deviation centred on the mean (Row et al., 2017; van Strien et al., 2012). We ranked models using the Bayesian information criterion (BIC) and report the top models within five BIC units (Richards, 2015). We used standardized beta coefficients to assess the direction of effect (positive or negative) of each landscape variable on genetic distance to infer how genetic distance is affected by landscape features and which variables have the greatest impact on neutral genetic variation.

Does divergent selection act across habitat types?

To identify loci associated with habitat type (pine vs. oak) while accounting for geography, we used a partial redundancy analysis (RDA; Forester et al., 2016) as our GEA test (hereafter referred to as GEA‐RDA), implemented in vegan (Oksanen et al., 2013). This constrained ordination approach models multivariate response data (in this case, genetic data) as a function of linear combinations of predictor variables by combining multivariate linear regression and PCA (Legendre & Legendre, 2012). RDA has been shown to be an effective method to detect weak, multilocus signatures of selection due to its low false positive and high true positive rates (Forester et al., 2018), making it an ideal test for GEAs in this system. We modelled our complete 3345 SNP data set (see Section 3) as a function of percent pine and percent oak within a 300‐m radius of each sampling location as habitat predictors of genetic divergence. We used the longitude of each individual's sampling location as a proxy for isolation‐by‐distance (see Section 3). Including longitude as a third “conditioned” matrix in the RDA allowed us to control for potential effects of isolation‐by‐distance. Longitude was weakly correlated with percent oak (|r| = .06) and percent pine (|r| = .35); therefore, we were able to include longitude in our model without confounding patterns between highly correlated variables. This allowed us to reduce potential noise driven by continuous population structure in our analysis of genetic divergence associated with habitat composition. We tested the significance of both the global model and model terms (percent pine and percent oak) using the anova.cca function in vegan with 10,000 permutations. Further, we identified outlier loci based on the “locus score,” which are the coordinates (loading) of each locus in the ordination space. We defined outliers as loci with loadings ±2.5 standard deviations from the mean on the two RDA axes (Forester et al., 2018). We then determined the predictor with which each outlier locus was most strongly associated based on the absolute value of the predictor's correlation coefficient for each SNP.

Does divergence in bill morphology have a genetic basis?

Langin et al. (2015) found that the magnitude of difference in bill length between pine and oak habitats in island scrub‐jays on Santa Cruz Island is similar to observed differences between pine/juniper and oak habitats in mainland populations of scrub‐jays (Peterson, 1993). The bill length difference among mainland scrub‐jays has also been shown to be adaptive for foraging on pine versus oak (Bardwell et al., 2001). Therefore, we chose to focus on bill length as our putatively adaptive trait of interest for identifying SNPs associated with bill morphology. Because island scrub‐jays show marked differences in bill morphology related to sex, age and season (Langin et al., 2015), it is noteworthy that one of the strengths of this data set is that all of our samples came from adult males that were measured by the same person (K. Langin) within the same 4‐month period, thereby reducing other sources of variation in the phenotypic data to increase the likelihood we would be able to detect loci associated with a complex morphological trait (Visscher et al., 2017). Bill length is positively associated with body size (Langin et al., 2015). Therefore, we performed a PCA on measurements of wing and tarsus length and extracted values from the first PC axis, which explained 68.9% of the variance, as an index of overall body size. We then used the residuals from a regression of bill length on this index of body size as a body size‐corrected measure of bill length in our analyses.

We used two different GWA tests to identify loci underlying variation in bill length using our imputed data set from beagle. We tested against the null hypothesis that bill morphology is entirely due to plasticity and therefore the marker effect size on bill length is zero. First, we used a Genome‐Wide Mixed Model Association algorithm (gemma; Zhou & Stephens, 2012) which fits a univariate linear mixed model using marker genotypes as a fixed effect in association tests with a single phenotype while correcting for population structure using a kinship matrix (random effect). Hereafter we refer to this approach as GWA‐GEMMA. By incorporating relatedness, we reduce the likelihood of spurious associations caused by neutral population structure (Sul et al., 2018). We applied a Benjamini–Hochberg correction for multiple testing to control for FDR and identified outlier SNPs based on a significance threshold of .05.

Univariate genome scan methods such as gemma have been shown to be biased towards large‐effect loci undergoing strong selection (Hoban et al., 2016; Wellenreuther & Hansson, 2016). Indeed, a Bayesian sparse linear mixed model showed that much of the significant genetic variation associated with bill length may be attributed to one SNP (see Supporting Information), which we would not necessarily expect, as bill morphology has been shown to be a complex polygenic trait in other systems (e.g., Bosse et al., 2017; Lundregan et al., 2018; Perrier et al., 2020). Our linear mixed model approach may therefore be biased towards these large‐effect SNPs while failing to detect the signal of smaller effect loci especially given our modest sample size. To address this limitation of gemma, we also used a partial RDA, which has been shown to be an effective method for identifying small‐effect loci (Forester et al., 2018), to model the effect of all 3345 loci on our body size‐corrected measures of bill length while controlling for population structure by using longitudinal coordinates of sampling localities as a proxy for isolation‐by‐distance (hereafter referred to as GWA‐RDA). Outlier loci were identified using the criterion of loading scores ±2.5 SD from the mean to identify SNPs associated with variation of bill length.

Visual inspection of linkage disequilibrium (LD) plots using LD calculated within a 100‐kb window in vcftools 0.1.17 with minimum and maximum alleles setting of 2 (Danecek et al., 2011) suggested high linkage between biallelic loci separated by ≤25 kb (Figure S2). We therefore input our zebra finch‐mapped SNPs in bedops version 2.4.39 (Neph et al., 2012), and output all genes within 25 kb of variant site coordinates in the annotated zebra finch genome. We then used the R packages mygene (Mark et al., 2020) and biomart (Durinck et al., 2005, 2009) to extract accession identifications and ontology information for each gene found within our query sequences using the available zebra finch Ensembl database. We compared all the loci flagged by our two GWA methods and searched annotation reports for GO terms related to bill morphology. We inferred those genes with these terms would best support our hypothesis of genetically based divergence of bill morphology.

RESULTS

Genotyping and filtering SNP matrices

We genotyped island scrub‐jays from oak habitat (n = 111), and western, central and eastern pine stands (n = 41; Figure 1) and obtained a total of 203,128,694 reads with an average of 1,336,373 reads per individual. Two individuals from the central oaks had low numbers of reads and were dropped from further analyses. Of the remaining reads, 91.8% were mapped to the Florida scrub‐jay reference genome with an average coverage of 16.8× per locus per individual following bwa‐mem alignment and genotyping in gstacks.

Initial processing of individuals with <24% missing data in stacks resulted in a matrix of 26,160 SNPs and 124 individuals. Two individuals were found to be closely related in the detect_duplicate_genomes step of radiator, so the individual with the least missing data in the pair was retained and the new matrix of 123 individuals and 25,815 SNPs was further filtered in radiator (Table S1). The number of SNPs and individuals kept after each filtering step are given in Table S1. This resulted in a final matrix of 3409 SNPs. pcadapt detected one putatively adaptive outlier locus, which was removed for a matrix of 3408 SNPs used for downstream neutral population genomic analysis (Section 3.2).

Because synteny is highly conserved in birds (Backström et al., 2008; Zhang et al., 2014), we were able to map 96.3% of our Florida scrub‐jay scaffolds to the zebra finch genome. This resulted in 3345 loci with chromosome‐level positional information for GEA (Section 3.3) and GWA (Section 3.4) analyses. pcadapt detected three putatively adaptive outlier loci within this imputed, mapped data set, which were removed for a matrix of 3342 SNPs for our PCAs (see below).

Population structure at neutral loci

We calculated the following population genomic parameters in stacks using the neutral SNP matrix of 3408 SNPs and 123 island scrub‐jays: observed heterozygosity (H O) = 0.256, expected heterozygosity (H E) = 0.272, inbreeding coefficient (F IS) = 0.066, and nucleotide diversity (π) across all sites (variant and fixed) = 0.00285. The effective population size (N e) estimate was 346.8 with a 95% confidence interval of 327–368 based on parametric bootstrapping. Mantel tests indicated significant (p = .001) spatial population structure for both measures of genetic distance (r PropShared = .235; r relatedness = .347) consistent with a pattern of isolation‐by‐distance primarily predicted by longitude according to our PCA results (Figures 2 and S3). We did not find evidence of discrete populations, and K = 1 was the best supported value of K by minimizing cross‐validation error in snmf (Alexander et al., 2009; Figure S4) and maximum log‐likelihood in structure (Pritchard et al., 2000; Table S2).

FIGURE 2

The island scrub‐jay population exhibited continuous spatial genetic structure across Santa Cruz Island. Coloured polygons represent woody vegetation where jays are found, with pine habitat outlined in black and the remaining shaded regions representing oak habitat. The colour gradient is a continuous representation of the predicted neutral genetic surface for the first principal component (PC) axis for sampled tree stands (grey = unsampled tree stands). Jays are expected to have a similar genetic composition if they are in areas with similar colours

MLPE analysis showed topographic distance and absolute difference in percent oak were consistently retained in top models and positively correlated with genome‐wide genetic distance of island scrub‐jays for both measures of genetic distance. Low vegetation density showed a negative correlation in models of the proportion of shared alleles, but was not retained in models of genetic relatedness. Absolute difference in percent pine was not retained in most models and showed a weak negative correlation with genetic relatedness (Table 1). None of our raw raster surfaces were highly correlated (|r| < .70) prior to running circuitscape. After calculating resistance surfaces, collinearity among environmental variables increased, due to circuitscape resistances being higher for individuals separated by greater geographic distances (McRae, 2006). Elevation was the only variable strongly correlated with other variables, with a VIF > 10, so it was removed from further analyses. Variables retained were topographic distance, low vegetation density, absolute difference in percent pine and absolute difference in percent oak (Figure S1).

TABLE 1

Maximum likelihood of population effects (MLPE) neutral landscape genomic results. Response variables were individual‐based genetic distances (proportion of shared alleles or relatedness). Fixed effects were topographic distance, low vegetation density, absolute difference in percent oak, and absolute difference in percent pine. Standardized beta coefficients (β) were used to assess the relative support for our a priori landscape resistance to gene flow hypotheses. The top supported models reported below are within 5 BIC units of the most supported model. Landscape resistance variables that had positive beta coefficients were considered to have a significant positive effect on genetic distance. Variables with negative beta interpreted here as nonsignificant

Genetic distance	Landscape feature	β	R 2	Delta BIC
Proportion shared alleles	Topographic distance	0.467	.273	0.0
	Low vegetation density	−0.179
	Absolute difference oak	0.047
	Topographic distance	0.478	.271	2.57
	Low vegetation density	−0.193
Relatedness	Topographic distance	0.428	.271	0.0
	Absolute difference oak	0.050
	Topographic distance	0.435	.272	1.57
	Absolute difference oak	0.066
	Absolute difference pine	−0.041

Evidence of divergent selection associated with habitat

The relative proportion of pine and oak within a 300‐m radius of each individual was significantly associated with genetic variation in island scrub‐jays based on RDA (p‐value percent pine = .001, p‐value percent oak = .001). A triplot of our two GEA‐RDA axes shows SNPs (dark grey points) and individuals (coloured circles) arranged in ordination space relative to their relationship with the predictor variables (black arrows, Figure 3a). Of the 168 outlier SNPs that were ±2.5 SD from the mean loading on RDA axis 1 (n = 93) and RDA axis 2 (n = 72; Figure 3b), we found 82 SNPs that were most strongly correlated with percent oak, and 86 SNPs that were most strongly correlated with percent pine (Figure S5).

FIGURE 3

Results of GEA and GWA of 123 island scrub‐jays. (a) Ordination plot of RDA showing the habitat‐linked (pine vs. oak) genetic divergence. Coloured points show where individual samples load for RDA axes 1 and 2 based on a partial RDA conditioned on the geographic location of each individual using 3345 SNPs (shown as dark grey points in the centre of the plot) as the response and relative proportion of pine and oak habitat within a 300‐m radius of sampling locality as the predictors (black vectors). The western, central and eastern panels represent individuals grouped by which pine stand they are geographically closest to. (b) Manhattan plot showing the absolute value of the SNP loadings on RDA axis 1 of the RDA‐GEA. (c, d) GWA results using univariate gemma and multivariate RDA, respectively, that identify SNPs associated with body size‐corrected measures of bill length. The horizontal dashed line in (c) corresponds to the threshold for statistical significance (p = .05). (d) Manhattan plot showing the absolute value of the SNP loadings on RDA axis 1 of the RDA‐GWA. The horizontal dashed lines in (b) and (d) show 2.5 SD from the mean absolute loading value for the RDA‐GEA and RDA‐GWA. Red highlighted points in (b) and (d) represent outlier SNPs associated with candidate genes (labels) related to bill morphology based on previous studies. Black and grey colours distinguish different chromosomes numbered according to the zebra finch nomenclature

Identification of loci underlying variation in bill morphology

Using GWA‐GEMMA, we detected seven SNPs that were significantly associated with bill length after correcting for multiple tests (FDR; Figure 3c), while the GWA‐RDA detected 91 additional SNPs (Figures 3d and S6). SNPs did not overlap between our two GWA analyses, though three of the GWA‐GEMMA SNPs were close to the GWA‐RDA outlier threshold, indicating some congruence of adaptive signatures across these different methods. Additionally, five SNPs detected by GWA‐RDA were also detected by our GEA‐RDA, three of which corresponded to candidate genes associated with bill morphology (COL14A1, PPP3CB, IGF1R; Figure 3; Table S3).

We identified 254 genes within 25 kb of our variant sites flagged by the two GWA analyses and GEA. Potential candidate genes found in our Gene Ontology term search of annotation reports (Table S4) included several genes previously identified in pathways associated with bill morphology (Abzhanov et al., 2004, 2006; Badyaev et al., 2008; Bosse et al., 2017; Lundregan et al., 2018; Mallarino et al., 2011) including bone morphogenetic protein (SMURF1; BMP15), ossification (PBX1), calmodulin (IQCB1; PPP3CB; HSPA2), Smad proteins (TRIM33), MAP kinase activity (“MAPK”: STK39; AKAP13; IGF1R), transforming growth factor beta (FN1), Notch (MFNG; TAFA4; TP63), and Wingless‐related integration (“Wnt”: ADGRA2; TRABD2B) signalling pathways (Figure 3b–d).

DISCUSSION

Divergence at a microgeographic scale

Evolutionary theory and empirical studies have long emphasized the importance of geographic isolation and reduced gene flow as a prerequisite for local adaptation to occur among populations. Nonetheless, local adaptation in the face of gene flow appears common even in taxa with a high dispersal potential (Richardson et al., 2014). For example, other avian systems have been shown to exhibit evolutionary responses in multiple morphological (Bertrand et al., 2016; Cornuault et al., 2015) and physiological traits (Gamboa et al., 2022; Mikles et al., 2020) within the potential dispersal radius of the study species. Understanding such fine‐scale adaptation requires multiple lines of evidence, including trait variation across environmental gradients and within the dispersal neighbourhood of the species in question, evidence for divergent selection across the environmental gradient, and observed phenotypic variation that must be at least partially genetically based.

Our study demonstrates that island scrub‐jays exhibit habitat‐linked microgeographic divergence. Neither snmf nor structure analyses revealed discrete population structure in island scrub‐jays (Table S2; Figure S4), suggesting that the spatial scale of Santa Cruz Island is within the dispersal capabilities of island scrub‐jays (Table 1; Figure 2). Our results indicate that despite this lack of structure, island scrub‐jays exhibited a genome‐wide pattern of isolation‐by‐distance (Table 1; Figure 2). We also observed relatively high genetic diversity comparable to estimates for mainland bird species (e.g., Mikles et al., 2020), and an N e estimate larger than found in other vertebrate species on Santa Cruz Island (Funk et al., 2016; Trumbo et al., 2021). Thus, localized isolation and genetic drift probably play a limited role in this system. Instead, we hypothesize that multiple, nonmutually exclusive processes including nonrandom gene flow and habitat choice (Edelaar & Bolnick, 2012), and/or divergent selection (Langin et al., 2015) are contributing to a pattern of microgeographic divergence in island scrub‐jays.

Localized dispersal is common across a wide range of taxa (Sexton et al., 2014) including Aphelocoma jays (Aguillon et al., 2017; Carmen, 1988; McCormack et al., 2008). Breeding island scrub‐jays maintain and defend year‐round territories (Caldwell et al., 2013; Collins & Corey, 1994; Curry & Delaney, 2020), and exhibit limited natal dispersal (Langin et al., 2015; Langin et al., 2017), but most individuals live as nonterritorial “floaters” for multiple years before a breeding territory becomes available (Curry & Delaney, 2020). Little is known about the movement patterns of jays during this transient period or about how dispersal decisions are affected by phenotypes or environmental context (Bolnick & Otto, 2013; Clobert et al., 2009). Similar to the Florida scrub‐jay, island scrub‐jays exhibit sex‐biased dispersal and may disperse up to 3 km from their natal territory (Langin et al., 2015). This distance is sufficient for the jays to move between adjacent oak and pine stands on Santa Cruz Island (Figure 1). Therefore, nonrandom dispersal via habitat selection might play a role in genetic differentiation even at this limited spatial scale (Edelaar & Bolnick, 2012; Edelaar et al., 2008; Nicolaus & Edelaar, 2018).

We found strong support for models containing multiple significant effects of landscape features on genetic distance in island scrub‐jays. Topographic distance was the strongest predictor of genetic distance, but we also observed that landscape features related to pine–oak ecotones affect genetic distance. The absolute difference in oak habitat at the scale of jay home‐range size was positively correlated with both our measures of genetic distance (Table 1). In other systems, studies have found support for both natal habitat preference (Camacho et al., 2016; Piper et al., 2013) and phenotype habitat‐matching (Benkman, 2016; Camacho et al., 2013; Garant et al., 2005). Both mechanisms have profound evolutionary consequences in the recognition of suitable habitat and settlement decisions (Berner & Thilbert‐Plant, 2015; Davis & Stamps, 2004; Edelaar et al., 2008) and are plausible hypotheses for reduced gene flow between habitat types in the island scrub‐jay. Our working hypothesis for the lack of an effect of absolute differences in pine habitat could be attributed to the fact that most jays included in this study were found primarily in oak habitat. Oak chaparral is more widespread on Santa Cruz Island compared to pine woodland (Junak, 1995), and most island scrub‐jays occur in oak habitat because of pine habitat being rare relative to oak habitat (Sillett et al., 2012). Consequently, we had a relatively small sample size of jays with substantial pine woodland within their home‐range (Figure 1). This may explain the low variation in pairwise differences in pine habitat presence among jay home‐ranges and could constrain our ability to detect a positive association between genetic distance and pine habitat (Forester et al., 2016; Landguth et al., 2012; Lotterhos & Whitlock, 2015). Our models also indicate that low vegetation density was not positively correlated with genetic distance (Table 1). Therefore, we do not have evidence to suggest these fully flight‐capable birds are incapable of dispersing across open areas despite their avoidance of marginal or barren habitats (Curry & Delaney, 2020).

We hypothesize that assortative mating could restrict gene flow across oak–pine ecotones and facilitate genome‐wide divergence between these habitat types (Servedio, 2016). Island scrub‐jays appear to mate assortatively based on bill morphology (Langin et al., 2015). While this pattern could be due to spatial autocorrelation in bill morphology across the island, there is some indication that active mate choice may play a role as well because female calls associated with pair‐bond formation and territory defense vary with bill morphology (Langin, Sillett, Morrison, et al., 2017; Podos 2001) . Future work investigating mate‐choice in island scrub‐jays will be needed to determine the mechanisms underlying apparent environmental barriers to gene flow in this system (Bradburd & Ralph, 2019; Manel et al., 2010).

We found support for divergent selection related to adaptation to pine versus oak habitats in island scrub‐jays. Individuals sampled in pine and oak habitat could clearly be separated based on their genotypes after controlling for isolation‐by‐distance (Figure 3a). These genetic signatures of selection are unlikely to be confounded by demographic processes due to the lack of discrete population structure in the island scrub‐jay (Figure S4; Cushman & Landguth, 2010; Lotterhos & Whitlock, 2015; de Villemereuil et al., 2014). Further, our association tests support the hypothesis that observed differences in bill morphology between island scrub‐jays living in pine and oak habitats have a genetic basis and are caused by many loci of small effect. We identified 96 unique SNPs out of our 3345 SNP data set that were strongly associated with variation in bill length using GWA‐GEMMA and GWA‐RDA, demonstrating at least a partial genetic basis to variation in bill morphology (Table S3). Notably, we found several candidate genes previously shown to be associated with variation in bill morphology, including those in the calmodulin pathway (IQCB1, PPP3CB, HSPA2; Abzhanov et al., 2006; Lundregan et al., 2018), and those involved in regulating BMP activity (SMURF1, FN1, IGF1R, AKAP13, MFNG, STK39, TRIM33; Abzhanov et al., 2004; Badyaev et al., 2008; Mallarino et al., 2011). One of the loci that was flagged by both our GEA‐RDA and GWA analyses was associated with the collagen trimmer, COL18A1 (Figure 3b,d), which is thought to be involved in palate development and is a significant marker for bill length (Bosse et al., 2017; Table S4). If adaptive developmental plasticity or simple differential wear between habitats were the only factors driving differences in bill morphology, we would not expect to see any evidence for loci associated with bill variation (Endler, 1986; Ghalambor et al., 2007; Kawecki & Ebert, 2004). In addition, we would not expect to find shared SNPs between our GEA‐RDA and GWA analyses, thereby linking selection associated with habitat to variation in bill morphology (Table S3).

Variations in bill morphology related to habitat gradients have been observed in several bird species and have been attributed to differences in local foraging substrates (Alonso et al., 2020; Badyaev et al., 2008; Benkman, 1993, 2003). Both bill depth and length have been shown to respond rapidly to selection due to food resource availability (Grant & Grant, 1993, 2002). The observed differences in bill length between island scrub‐jays living in pine and oak habitats mirror patterns seen in mainland Aphelocoma jays, where bill shape varies adaptively with local resources, but on a spatial scale of a few kilometers instead of hundreds of kilometres (Bardwell et al., 2001; Langin et al., 2015; McCormack & Smith, 2008; Peterson, 1993). Although we have not quantified foraging behaviour or other ways the bill is used in island scrub‐jays within pine and oak habitats, prior research of mainland jays shows that bill morphology significantly affects foraging efficiency on different substrates (Bardwell et al., 2001), which may influence fitness (Pyke, 1984). Additional research quantifying the feeding performance on different diets relative to bill morphology will be needed to understand the potential fitness consequences of phenotype–habitat mismatch.

While not the focus of this study, our GEA analysis also flagged multiple SNPs correlated with habitat type that were associated with genes without a direct relationship to bill morphology (Tables S3 and S4). Previous research has found genetic differentiation across fine spatial scales in several avian species related to physiological and reproductive adaptive divergence between habitats (e.g., Charmantier et al., 2016; García‐Navas et al., 2014; Mikles et al., 2020; Perrier et al., 2020; Senar et al., 2006). We hypothesize that the signature of divergent selection at these loci could be attributed to environmental differences between habitats driven by the climate gradient induced by the California Current: Santa Cruz Island is characterized by more arid conditions in the east and a cooler climate in the west (Fischer et al., 2009; Gamboa et al., 2022; Morrison et al., 2011). Given that birds exhibit strong physiological responses to extreme heat (Mckechnie & Wolf, 2009), the island's pronounced east–west climate gradient could lead to divergence in physiological or behavioural traits in the island scrub‐jay. Hence, our observed pattern of selection across habitats may involve adaptive differences in traits other than bill length.

Conservation implications for a range‐restricted species

Adaptability and resilience of wild populations, particularly those with restricted ranges, are threatened by rapid climate change, emphasizing the need to understand the evolutionary processes underlying genetic and phenotypic microgeographic variation (Funk et al., 2019; Hohenlohe et al., 2021). The island scrub‐jay is a species of conservation concern because it currently occurs on a single, small island (Morrison et al., 2011), is a key seed disperser (Pesendorfer et al., 2016), and is vulnerable to catastrophic population declines related to climate change and disease (Bakker et al., 2020). Given that the adaptive divergence in island scrub‐jays living in pine and oak habitats is genetically based, conservation planning should consider the implications of such microgeographic variation in management decisions (Langin et al., 2015). For example, instead of only using abundance as a metric of population viability, managers could incorporate information on adaptive variation driven by habitat heterogeneity to effectively conserve the adaptive capacity of individuals across their range; such capacity will probably only increase in importance as climates rapidly change (Funk et al., 2019, 2012; Hohenlohe et al., 2021).

CONCLUSION

The island scrub‐jay exhibits a remarkable repeated pattern of habitat‐linked differences in bill morphology between adjacent pine and oak stands hundreds of metres apart. Prior evidence for a genetic basis to this pattern was based solely on pedigree‐based heritability estimates that were potentially confounded by spatial autocorrelation in bill morphology (Langin et al., 2015). In this study, we characterized fine‐scale genetic variation using genome‐wide SNP data to test for genomic signatures underlying this pattern of microgeographic divergence. Neutral landscape genomic analyses showed a significant pattern of isolation‐by‐distance and increased genetic distance associated with habitat variation consistent with divergence within the dispersal capabilities of island scrub‐jays. In addition, we found signatures of divergent selection at the genomic level by identifying multiple SNPs from across the genome that were strongly associated with habitat composition and bill length across the pine–oak ecotone, and evidence that bill morphology has a genetic basis. Collectively, our results support the hypothesis that habitat‐linked divergence in bill morphology in island scrub‐jays represents microgeographic adaptation within a single population. The island scrub‐jay system provides a robust framework to further test for microgeographic adaptation, and in so doing, also to inform the conservation of one of North America's most range‐restricted bird species.

CONFLICT OF INTEREST

The authors declare they have no conflict of interests.

AUTHOR CONTRIBUTION

W.C.F., C.K.G., T.S.S., and K.M.L. designed the study. K.M.L collected field data with support from S.A.M. and T.S.S.. P.E.S. designed and conducted laboratory work and developed genomic libraries. R.G.C. analysed data with guidance and assistance from W.C.F., B.R.F., P.E.S., D.R.T and N.C. The article was written by R.G.C. with input from all authors.

OPEN RESEARCH BADGES

This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.8sf7m0cpq, and https://github.com/RGCheek/ISSJ_Genomics_150_RAD.

Fig S1‐S6

Click here for additional data file.

Table S1‐S4

Click here for additional data file.