Environmental Drivers of Genetic Structure and Local Adaptation in a Marine Foundation Species

TL;DR

This study explores how environmental factors influence genetic diversity and adaptation in a kelp species across a large geographic range, revealing patterns of local adaptation and genetic differentiation.

Contribution

The study provides new insights into the genetic structure and local adaptation of Eisenia arborea in response to environmental gradients and climate change.

Findings

Strong genetic differentiation was found between northern and southern populations of Eisenia arborea.

Environmental association analyses identified SNPs correlated with sea surface temperature and depth, indicating local adaptation.

Southern populations showed low genetic diversity and high inbreeding, suggesting vulnerability to climate change.

Abstract

Predicting how species will respond to global change requires understanding how environmental drivers shape both neutral and adaptive genetic variation across space. The kelp
Eisenia arborea
is a thermally tolerant foundation species spanning more than 3000 km of coastline and a broad latitudinal temperature gradient in the Northeast Pacific, yet how environmental and demographic processes influence genomic and population structure remain unclear. We used genome‐wide ddRAD sequencing to investigate patterns of genetic diversity, connectivity and local adaptation in
E. arborea
across two depths and ~2700 km of coastline. We detected strong genetic differentiation between northern (British Columbia, Canada) and southern (Baja California, Mexico) populations, indicating limited gene flow across the species' broad range. Southern populations also had the lowest genetic diversity and highest inbreeding, broadly consistent with expectations for populations occupying environmentally marginal portions of a species' range. However, the two southernmost populations (~200 km apart) were highly similar and well connected, whereas mid‐range sites were more differentiated, indicating that the geographic range edge population was not genetically isolated as is often hypothesised. Environmental association analyses identified SNPs correlated with both sea surface temperature and depth, revealing signals of local adaptation to broad climatic gradients and fine‐scale habitat variation. The combination of high inbreeding, restricted connectivity and local adaptation highlights both the vulnerability and potential conservation value of distinct genetic units, especially warm‐adapted southern populations, for maintaining the resilience of these Eisenia forests under ocean warming.

We used genome‐wide SNP data to investigate how environmental gradients shape genetic structure, connectivity and local adaptation in the kelp
Eisenia arborea
across ~2700 km of coastline and contrasting depth environments. We found strong latitudinal genetic structure, evidence for local adaptation to temperature and depth despite high connectivity within locations, and demographic signatures consistent with marginality at the warm range edge. Together, these results highlight how neutral and adaptive processes interact across spatial scales in a marine foundation species, with implications for conservation under ongoing climate change.

Figures

Genetic structure of Eisenia arborea populations across the west coast of North America. (a) Principal component analysis (PCA) showing clustering of individuals from each location. (b) sNMF structure plot showing population structure (k = 4) of all individuals. (c) A map showing the location of the six sampling locations relative to the known distribution of E. arborea . Each black point is a field observation, compiled from https://macroalgae.org while each yellow point indicates a sampling site.

Population connectivity and differentiation of Eisenia arborea . (a) Relative migration (based on Jost's D method) between sites (those with less than 10 samples are omitted). Each circle represents a site and the inferred relative migration is shown by the arrows. Only values greater than 0.05 are shown; all other values were between 0.010 and 0.049. (b) Pairwise F ST values calculated between sites coloured by their values. High F ST indicates high differentiation between sites.

Genetic diversity statistics and effective population size estimates across samples. (a) Inbreeding coefficient, (b) observed heterozygosity and (c) expected heterozygosity for each sample of E. arborea .

Demographic reconstructions using SFS implemented in StairwayPlot2. (a–g) Effective population size (N e) through recent geological time for each population with sufficient sampling (min. 16 haplotypes following downscaling in easySFS). (h) Effective population size (log‐transformed) at time zero (i.e., present) for each population as inferred from StairwayPlot2. Grey shading indicates the 95% confidence interval of inferred population size.

Ordinations of neutral versus potentially adaptive SNPs. (a) PCA of neutral SNPs (those not identified as outliers of environmentally‐associated SNPs). (b) PCA of loci significantly associated with an environmental variable (EAA) or identified as an outlier from Bayescan.

Loci potentially under selection across environments. (a) Examples of loci significantly associated with maximum SST and reflecting underlying latitudinal structure. (b) Examples of depth‐associated SNPs.

Tables

TABLE 1. Sample sizes and private alleles associated with each site and population.

Population	Cluster	Latitude	Depth	Sample size	N p	N p (loc.)	N p (cluster)
Bahía Magdalena	1	24.5	Shallow (1 m)	N = 21	181 (0.02)	460 (0.05)	2130 (0.23)
	1		Deep (10 m)	N = 5	35 (< 0.01)
Bahía Asunción	1	27.1	Shallow (1 m)	N = 20	438 (0.05)	1069 (0.11)
	1		Deep (10 m)	N = 8	84 (< 0.01)
San Quintin	2	30.6	Shallow (3 m)	N = 10	58 (< 0.01)	477 (0.05)	477 (0.05)
	2		Deep (10 m)	N = 13	81 (< 0.01)
Ensenada	3	31.8	Shallow (2 m)	N = 18	435 (0.05)	1724 (0.19)	1724 (0.19)
	3		Deep (12 m)	N = 16	267 (0.03)
Taylor Islet	4	48.8	Shallow (1 m)	N = 7	54 (< 0.01)	54 (< 0.01)	1323 (0.14)
Danver's Point	4	48.9	Shallow (1 m)	N = 14	145 (0.02)	1262 (0.14)
	4		Deep (8 m)	N = 6	74 (< 0.01)

Full text

Introduction

1

Understanding how organisms adapt to environmental conditions is fundamental for predicting their biological responses to global change (Atkins and Travis 2010; DeMarche et al. 2019). While many studies have explored how broad‐scale climatic gradients shape genetic structure (e.g., Coleman et al. 2011; Fabian et al. 2012; Zhang et al. 2019; Vranken et al. 2021), we lack a clear understanding of how fine‐scale environmental variation interacts with large‐scale drivers to shape genetic diversity and local adaptation (López‐Goldar and Agrawal 2021; Marrot et al. 2022). This remains a key gap because such interactions can either amplify or buffer the effects of climate change, influencing both short‐term persistence and long‐term evolutionary potential (Cortés 2017; DeMarche et al. 2019; Starko, van der Mheen, et al. 2024). Fine‐scale gradients in key environmental traits, such as light or resource availability, can act as critical drivers of natural selection that co‐occur with broader‐scale gradients, such as latitudinal temperature variation. Moreover, latitudinal and local‐scale drivers can interact to lead to complex effects on species performance (Giraldo‐Ospina et al. 2020; Baum et al. 2023; Starko et al. 2023, 2025), suggesting that the effects of natural selection in the face of multiple drivers may be complex and challenging to predict. This is especially important in coastal marine environments, where steep environmental gradients can occur over short spatial scales (Harley et al. 2006; Starko, van der Mheen, et al. 2024) and potentially drive adaptive divergence even in the presence of high gene flow (e.g., Coyer 2004; Minne 2024).

At broader spatial scales, the central‐marginal hypothesis (CMH) provides a useful framework for interpreting variation in genetic diversity across a species' range (see synthesis by Eckert et al. 2008). CMH predicts that populations towards the geographic edges of a species' distribution may exhibit reduced genetic diversity, smaller effective population sizes, and greater genetic differentiation due to isolation and demographic instability that arises as suitable habitat becomes rarer towards the range edge of a species (Eckert et al. 2008). However, a key assumption of CMH is that geographic range edges correspond to environmentally marginal conditions and that habitat is more optimal and/or more abundant at a species' range center, an assumption that may not always hold in heterogeneous environments (Micheletti and Storfer 2015; Pironon et al. 2017; Macrì et al. 2021). In marine systems, in particular, environmentally marginal or stressful conditions can occur in a patchy manner across a species' range. Thus, patterns of genetic diversity and differentiation may reflect a combination of geographic range position and local or regional environmental contexts. Although population genetic patterns in some marine species are consistent with CMH predictions (e.g., Alberto et al. 2006; Liggins et al. 2014; Vranken et al. 2021), others have reported spatial patterns that deviate from expectations (e.g., Kennedy et al. 2020; Gierke et al. 2023; Ha et al. 2024; Edgeloe et al. 2025), highlighting variability across taxa and different environments. Understanding how range dynamics influence genomic patterns across space requires examining not only neutral genetic patterns shaped by demographic history, but also how spatial variation in selection driven by environmental gradients interacts with those demographic processes. Considering both demographic (neutral) and adaptive (selective) perspectives provides a more complete picture of the forces shaping population structure and may help explain inconsistencies in how range position influences genomic patterns among marine species.

Kelp forests, formed by large brown seaweeds primarily from the order Laminariales, are among the ocean ecosystems most threatened by climate change (Pörtner et al. 2019; Wernberg et al. 2024). As kelp populations fail to adapt to rapidly increasing stressors, such as marine heatwaves which are becoming longer, more intense and more frequent with climate change (Oliver et al. 2018, 2019), large areas of kelp forest are being rapidly lost (Filbee‐Dexter et al. 2022; Wernberg et al. 2024). This is particularly problematic given that only limited restoration success has been documented for these habitats (Wood et al. 2024). Moreover, climate change‐amplified marine heatwaves are already restructuring the genetic make‐up of habitat‐forming seaweeds and, in some cases, eroding their genetic diversity (Coleman, Minne, et al. 2020; Gurgel et al. 2020), potentially influencing the adaptive capacity of species. While there have been growing efforts to understand spatial patterns of kelp genetic diversity and structure (e.g., Bringloe et al. 2021; Vranken et al. 2021; Bemmels et al. 2025) and how these patterns relate to broad‐scale environmental drivers (e.g., Coleman et al. 2011; Vranken et al. 2021; Minne et al. 2025), we have a very limited understanding of how fine‐scale environmental variation interacts with broader climatic drivers (e.g., temperature variation) to influence genetic diversity and structure of natural kelp populations. For example, it remains unclear whether selection imposed by local‐scale environmental drivers can meaningfully influence the genetic characteristics of kelp populations even in the face of high levels of gene flow and stronger selective drivers acting on broader scales (e.g., latitudinal climate variation). Given that the impacts of warming on marine foundation species can vary substantially across different genotypes (Burgess et al. 2021; Starko et al. 2023; Starko, van der Mheen, et al. 2024), this remains a critical knowledge gap when trying to predict the capacity of species to respond to and adapt to climate change and heatwaves.

Eisenia arborea (J. E. Areschoug 1876) is an ecologically significant stipitate and palm‐like kelp species that is endemic to the Pacific coast of North America. The geographical distribution of E. arborea extends from Baja California Sur (Mexico) to Haida Gwaii (Canada) (Druehl 1970), although it exhibits an atypical disjunct distribution where it is common in Baja California and Southern California, absent between Central California and British Columbia (Druehl 1970), and then moderately widespread but patchy in British Columbia (Watson et al. 2021). E. arborea has been observed to thrive in a variety of marine habitats, ranging from intertidal rocky substrates to depths exceeding 30 m. This species can form canopy‐like structures, both as dominant habitat‐formers and as patches in the understory of giant kelp ( Macrocystis pyrifera ) forests (Schiel and Foster 2015), thus it provides vital habitats for diverse marine communities and commercially valuable species (Serviere‐Zaragoza et al. 1998; Sievers et al. 2016).

E. arborea is distributed more equatorward (i.e., south) on the mainland than any other kelp species in North America (Druehl 1970; Avila‐Pedroche et al. 2008; Hernandez‐Carmona et al. 2011), and therefore it is likely the most thermally tolerant kelp in the Northeast Pacific. In the southernmost regions of the Baja California peninsula, E. arborea populations have demonstrated the capacity to cope with temperatures considered very warm for kelp species (e.g., 23°C–25°C or more) and low nitrate availability (often less than 1 μM) (Matson and Edwards 2007; Hernandez‐Carmona et al. 2011; Zertuche‐González et al. 2022). In contrast to other kelps, such as Macrocystis pyrifera , E. arborea sporophytes have been shown to exhibit resilience to anomalous warming events associated with El Niño Southern Oscillation (ENSO) (Ladah et al. 1999; Edwards and Hernández‐Carmona 2005) and other drivers of marine heatwaves (Arafeh‐Dalmau et al. 2019). Moreover, experimental work has demonstrated high thermal tolerance of the early life stage of E. arborea (Matson and Edwards 2007; Muth et al. 2019; Almeida‐Saá et al. 2025). This may explain not only why populations appear resilient in southern parts of the Baja California peninsula, but also why they appear to be expanding in parts of British Columbia (Watson et al. 2021), as warming threatens other kelp species (Starko et al. 2022, 2025; Starko, Timmer, et al. 2024; Mora‐Soto et al. 2024).

The perennial and stress‐tolerant nature of E. arborea , along with its favourable chemical composition for several industrial purposes (e.g., fucoidan and alginate production), have fueled commercial harvesting and mariculture of this species (Zertuche‐González et al. 2014, 2022), sustaining socio‐economic benefits in Baja California (Landa‐Cansigno et al. 2017; Múzquiz de la Garza et al. 2019). Despite the economic and ecological importance of E. arborea , its genetic diversity and population structure remain poorly understood, limiting our understanding of the adaptive capacity of populations or how best to regulate restoration or aquaculture activities (Coleman, Wood, et al. 2020). However, prior genetic and phenotypic studies have found evidence of local adaptation in this species to environmental conditions, such as nutrient availability and hydrodynamic forces, over short distances (Coyer 2004; Matson and Edwards 2006; Parada et al. 2012). Consequently, E. arborea may offer an interesting opportunity to understand local adaptation in a habitat‐forming kelp and how adaptation to fine‐scale environmental variation might vary across broad latitudinal gradients. For example, kelps growing at different depths will experience different light levels and may experience differences in thermal variability with higher fluctuations towards the surface (Bongaerts et al. 2010; Giraldo‐Ospina et al. 2020; Starko, van der Mheen, et al. 2024). This depth‐driven environmental variability is likely to act as a source of selection, potentially altering the allele frequencies of key genes associated with light acquisition or thermal tolerance (Minne 2024). However, resolving whether local‐scale environmental drivers can meaningfully influence the genetic characteristics of kelp populations is essential, because it determines the spatial scale at which management and conservation interactions will be most effective (Starko, van der Mheen, et al. 2024).

In this study, we assessed the genetic diversity, connectivity and population structure of both deep (8–12 m) and shallow (1–3 m) Eisenia arborea stands across a broad latitudinal range spanning more than 800 km in Mexico as well as two locations in Barkley Sound, British Columbia, ~1925 km north of the highest latitude Mexican sampling site. Sampling along both depth and latitudinal gradients allowed us to evaluate the influence of large‐scale climatic gradients and small‐scale environmental heterogeneity on the genetic structure of this foundation species. Specifically, we leveraged this genomic dataset to test whether genetic structure and diversity were consistent with the CMH (Eckert et al. 2008); specifically whether southern range edge populations have reduced genetic diversity, elevated inbreeding, increased genetic differentiation and historically reduced effective population sizes, a pattern seen in some past studies on large habitat‐forming seaweeds (e.g., Clark et al. 2020; Vranken et al. 2021), but not others (e.g., Gierke et al. 2023; Edgeloe et al. 2025). We also tested the hypothesis that local environmental gradients (in this case depth) can offer sufficient selection to alter allele frequencies despite high gene flow and the presence of broader drivers of selection, such as latitudinal climate variation. Importantly, depth serves as a proxy for fine‐scale environmental heterogeneity that integrates multiple covarying factors including light availability, temperature variability, nutrient flux and hydrodynamic exposure to waves and currents. Importantly, the relative importance and covariance of these factors are expected to vary among locations, such that selection associated with depth may be context dependent rather than uniform across the species' range. Our findings provide important insights into the factors that structure the genetic diversity of marine species with specific implications for understanding the resilience and adaptive capacity of E. arborea populations in the face of ongoing environmental change.

Methods

2

Field Sampling

2.1

We sampled kelp from locations (n = 6) spanning the latitudinal distribution of E. arborea (Figure 1). At all but one location, we collected from two depths (see Table 1 for collection information), leading to 11 sampled sites. Samples of laminae were collected and kept cool (usually in a cooler with ice) until they could be processed, at which point they were cut into thin pieces using a razor blade or scalpel and added to Qiagen AE buffer for storage until further processing. All sampling was conducted via snorkel or SCUBA. While we aimed for 15–30 individuals per site, several samples were not viable for genomic‐scale sequencing (see below) due to issues with DNA quality, causing final sample sizes to vary from 5 to 21 individuals across sites (Table 1).

Genetic structure of Eisenia arborea populations across the west coast of North America. (a) Principal component analysis (PCA) showing clustering of individuals from each location. (b) sNMF structure plot showing population structure (k = 4) of all individuals. (c) A map showing the location of the six sampling locations relative to the known distribution of E. arborea . Each black point is a field observation, compiled from https://macroalgae.org while each yellow point indicates a sampling site.

DNA Extraction, Library Preparation and Sequencing

2.2

Total genomic DNA was extracted using the Qiagen DNeasy Plant Mini Kit, following a modified protocol (Vranken et al. 2021). Specifically, the incubation period was extended to 48 h, and an additional wash with 95% ethanol was conducted before DNA elution from the spin column. A DNA clean‐up was then performed prior to library preparation using the Qiagen PowerClean kit. DNA concentration was measured with a Qubit fluorometer (Invitrogen), and quality was verified using either a LabChip GX Touch 24 (PerkinElmer) or standard gel electrophoresis. A total of 200 ng of DNA per sample was used to create ddRAD (double digest restriction associated DNA) libraries according to Severn‐Ellis et al. (2020), with some modifications. After digestion with PstI and NlaIII restriction enzymes (New England Biolabs), unique barcoded adaptors were attached to the DNA fragments of each sample. Double size selection was conducted with SPRI beads (AMPure XP, Beckman Coulter) to retain fragments within the 250–800 bp range. Libraries were amplified by polymerase chain reaction (PCR) and purified with SPRI beads. Equimolar concentrations of 48 individual libraries were pooled to create pooled libraries, with up to six pooled libraries sequenced on an Illumina NovaSeq X lane as part of a larger project.

SNP Calling and Filtering

2.3

Paired‐end reads were demultiplexed using the process_radtags function in STACKS (2.64) (Catchen et al. 2013), applying quality filtering (‐q, ‐c), barcode recovery (‐r) and RADtag verification (‐‐renz_1 pstI; ‐‐renz2 NIaIII) following Vranken et al. (2021). Single nucleotide polymorphisms (SNPs) were then identified de novo in STACKS. Adapter trimming was performed using TRIMMOMATIC with the parameters ILLUMINACLIP:TruSeq3‐PE.fa:2:30:10 before all reads were trimmed to 140 bp with BBMAP (BBMap—Bushnell B.—sourceforge.net/projects/bbmap/T). Sequence read quality was assessed with FASTQC (Andrews 2010) and MULTIQC (Ewels et al. 2016). The denovo_map pipeline was used for assembly and SNP calling, with a minimum stack distance of three nucleotides (‐m) and a maximum distance of three nucleotides allowed between stacks within a locus (‐M). Up to three mismatches were permitted among orthologous loci across individuals in the catalogue. To filter low‐quality SNPs, VCFTOOLS version 0.1.1671 was used to remove all indels and sites with more than two alleles. We set a minimum coverage depth of 5, a maximum coverage at twice the mean sequencing depth, and excluded sites with over 10% missing data or a minor allele frequency below 0.02. We also removed all individuals with missing data from more than 10% of loci. We used technical replicates (n = 8 pairs) to calculate the sequencing error rate. For 75% (6/8) of samples, error rates were below 1% while two replicates had error rates of 1.3% and 3.1%, respectively (Figure S1). The two replicate pairs with higher error rates were those with the most missing data, suggesting that the higher error rates reflect lower quality or coverage of sequence reads.

Data Analysis

2.4

Most other data manipulation and analyses were carried out in R 4.3.2 (R Core Team 2023). Data manipulation involved the packages vcfR (Knaus and Grünwald 2017), tidyverse (Wickham et al. 2019), dartR (Gruber et al. 2018), SNPRelate (Zheng and Zheng 2013) and ape (Paradis and Schliep 2019). To identify population structure, we conducted principal component analysis (PCA) in vegan (Oksanen et al. 2010) and structure analysis in sNMF using LEA (Frichot and François 2015), selecting the number of clusters (K) to equal that which leads to the lowest cross‐entropy or to the value beyond which reductions in cross‐entropy are minimal (i.e., the ‘elbow’ method). We quantified genetic differentiation among sites using both pairwise F ST and Jost's D. While F ST remains widely used, it can underestimate differentiation when within‐population heterozygosity is high. Jost's D provides a complementary measure of allelic differentiation that is less dependent on within‐population diversity and has been shown to perform well for SNP datasets (despite originally being developed for multiallelic markers) when interpreted as a relative measure of differentiation rather than a direct estimator of migration (Zimmerman et al. 2020). We therefore interpret both metrics as complementary indicators of relative genetic differentiation. We calculated pairwise F ST and various genetic diversity metrics in hierfstat (Goudet 2005) including significance testing using permutation tests (n = 10,000 permutations), in which individuals were randomly reassigned among populations. p‐values were adjusted for multiple comparisons using the Bonferroni method. We then used the SNP data to infer relative migration using the divMigrate function from the diveRsity package (Keenan et al. 2013) implementing Jost's D. For relative migration, we only included sites with 10 or more individuals to avoid issues associated with very low sample sizes. We tested for isolation by distance using a Mantel test implemented in adegenet comparing geographic distance obtained using the geosphere package (Hijmans et al. 2017) and F ST/(1 − F ST) values between paired sites.

To identify SNPs that were putatively under selection, we ran two environmental association analyses (EAA) as well as an F ST‐based outlier detection approach. Both EAA analyses and outlier detection were conducted using all 9315 SNPs. Our goal with EAA was to identify SNPs that were associated with both the broad‐scale temperature gradient present latitudinally and with local‐scale depth gradients. Thus, we used only annual maximum sea surface temperature (SST max) from Bio‐Oracle 3 (Assis et al. 2024) and field‐measured depth (two levels = deep, shallow). We chose to only include maximum SST at the exclusion of other temperature measurements because maximum temperature is considered to be a major driver of kelp distributions (Wernberg et al. 2024) and due to high covariance between maximum, minimum and average sea surface temperature data. For EAA analyses, we used both redundancy analysis (RDA) and latent factor mixed models (LFMM) implemented in R 4.3.2 (R Core Team, 2021), recognising that these approaches are distinct statistical methods (multivariate vs. univariate) with different sensitivities. Rather than restricting our analyses to only the overlapping loci between these methods (as is sometimes done in landscape genomics studies), we retained all loci identified as significant by either approach. This decision was based on growing evidence that different EAA methods may identify complementary subsets of loci under selection, particularly when the strength and consistency of selection varies across space (Capblancq et al. 2018; Forester et al. 2018; Capblancq and Forester 2021), which we expected to be the case with our dataset. We also conducted outlier analysis in BayeScan (Foll 2012) to identify putative loci under selection using an F ST based approach.

We conducted reconstructions of historical population size changes through geological time using site frequency spectra implemented in both StairwayPlot2 (Liu and Fu 2020) and EPOS (Lynch et al. 2020). In order to conduct these analyses only on neutral loci, we removed all SNPs that were identified with any EAA or outlier method. We also randomly selected one SNP per ddRADtag (i.e., thinning) to avoid non‐independence of closely linked loci. For calculation of SFS, we used a downscaling method to account for missing data using the easySFS script (https://github.com/isaacovercast/easySFS) in Python 3.12.3 (Python Software Foundation). We only included sites to which we could reasonably downscale to at least 16 haplotypes (i.e., 8 diploid individuals) to avoid error associated with low sample sizes. Downscaled sample sizes (number of haplotyes) were as follow: Magdalena Bay Shallow: 26; Bahia Asuncion Shallow: 30; San Quintin Shallow: 16; San Quintin Deep: 16, Ensenada Shallow: 28; Ensenada Deep: 22; Danver's Point Shallow: 18. We note that even with sufficient sample sizes, demographic analyses should be interpreted with caution when using fewer than 10^6^ SNPs (Lapierre et al. 2017).

Results

3

After filtering, our resulting dataset contained 9315 single nucleotide polymorphisms (SNPs) across 138 individuals (see Table 1). Across all sites, we identified four distinct clusters using principal component analysis (Figure 1a,b). Perhaps not surprisingly, deep and shallow samples from the same location clustered together in every case, as did the two locations in British Columbia, which were only ~10 km apart, indicating that nearby sites shared high genetic similarity. However, kelp from the two southernmost locations (Bahía Magdalena and Bahía Asunción) was also found to cluster together despite the nearly 200 km distance between these locations. Based on cross‐entropy criteria from the sNMF structure analysis, the optimal number of ancestral lineages was four (K = 4). However, we also plotted K = 5 (reflecting a separation of Bahía Magdalena and Bahía Asunción locations; Figure S2) which had a slightly lower cross‐entropy but only marginally lower than K = 4 (Figure S3) and so it was considered suboptimal.

Population connectivity and differentiation metrics also supported these patterns of population structure. Pairwise differentiation based on F ST values was statistically significant for all population comparisons (permutation tests, p < 0.01), indicating subtle yet consistent genetic differentiation across the species' range (Figure 2b). Despite universal statistical significance, the magnitude of differentiation varied strongly among site‐level comparisons (0.01–0.94), with the lowest values observed between depth strata within locations and between the two southernmost populations, and highest differentiation among populations at opposite ends of the geographic range of the species (e.g., Magdalena versus Barkley Sound). This pattern led to a strong, significant pattern of isolation‐by‐distance (Mantel test: r = 0.8761, p < 0.001; Figure S4). While there was high connectivity inferred from Jost's D (Figure 2a) across different depths at the same location, most of the sampling locations were meaningfully differentiated based on both F ST and Jost's D (Figure 2). The exceptions here are that, again, the two British Columbia locations were highly similar according to F ST and that the two southern‐most locations in Mexico showed high connectivity, based on Jost's D, and high similarity based on F ST.

Population connectivity and differentiation of Eisenia arborea . (a) Relative migration (based on Jost's D method) between sites (those with less than 10 samples are omitted). Each circle represents a site and the inferred relative migration is shown by the arrows. Only values greater than 0.05 are shown; all other values were between 0.010 and 0.049. (b) Pairwise F ST values calculated between sites coloured by their values. High F ST indicates high differentiation between sites.

Genetic diversity was lowest and inbreeding was highest towards the southern range edge of E. arborea . Heterozygosity (H_E_ and H_O_) was consistently lower at the two warmest locations compared to all other locations and this was consistent across depths (Figure 3b,c). Moreover, inbreeding coefficients (F IS) were high across all sites, ranging from 0.23 to 0.70 but tended to be highest towards the warm range edge (Figure 3a; Table S1). In particular, the shallow site from Bahía Magdalena and the deep site from Bahía Asunción had the highest inbreeding coefficients of any sites that we sampled. Sites in British Columbia, however, especially the deep Danver's site, also had notably high F IS values.

Genetic diversity statistics and effective population size estimates across samples. (a) Inbreeding coefficient, (b) observed heterozygosity and (c) expected heterozygosity for each sample of E. arborea .

Analysis of site frequency spectra (SFS) suggested that most populations have experienced declines in effective population size over extended evolutionary timescales (Figure 4, Figure S5). Given the absence of species‐specific estimates of mutation rate, temporal estimates inferred from SFS‐based analyses are highly uncertain and should be interpreted qualitatively. Nonetheless, across populations, SFS patterns were generally consistent with a history of larger ancestral population sizes followed by declines towards the present. Shallow populations from Magdalena and Danver's Point exhibited lower inferred historical Ne than most mid‐range sites, whereas populations from San Quintín showed comparatively more stable demographic trajectories in recent portions of the reconstruction (Figure 4).

Demographic reconstructions using SFS implemented in StairwayPlot2. (a–g) Effective population size (N e) through recent geological time for each population with sufficient sampling (min. 16 haplotypes following downscaling in easySFS). (h) Effective population size (log‐transformed) at time zero (i.e., present) for each population as inferred from StairwayPlot2. Grey shading indicates the 95% confidence interval of inferred population size.

EAA analyses identified numerous loci associated with both SSTmax and depth (n = 1368 total; Figure S6) as well as several outlier loci according to BayeScan. However, very few potentially depth‐associated loci were identified by more than one detection method. We identified 137 environmentally‐associated SNPs with RDA (depth = 93, SST = 44), 1231 SNPs with LFMM (depth = 28, SST = 1203), and 31 outliers with BayeScan. Ordination of these loci (Figure 5b) reveals similar patterns to those from all SNPs (Figure 1a) or the neutral SNPs (Figure 5a) with some key differences. While neutral loci reflect only four clear clusters, the EAA and outlier SNPs (hereafter ‘putatively adaptive loci’) show clear separation of each location with six evident clusters. Moreover, while Ensenada and San Quintin are quite distinct in the neutral SNPs, they are similar with respect to the putatively adaptive SNPs, potentially reflecting selection from similar environments due to the relatively close proximity of these locations.

Ordinations of neutral versus potentially adaptive SNPs. (a) PCA of neutral SNPs (those not identified as outliers of environmentally‐associated SNPs). (b) PCA of loci significantly associated with an environmental variable (EAA) or identified as an outlier from Bayescan.

While we identified a much stronger effect of temperature, several SNPs were significantly depth‐associated despite the high levels of gene flow across depths at the same location. Temperature‐associated loci were generally quite consistent, reflecting a transition from one allele to another along the broad latitudinal gradient. However, this transition happened at a range of locations along the north–south gradient depending on the SNP (Figure 6a). In contrast, loci associated with depth reflected a broader range of patterns. Many of these loci (e.g., locus 5844, locus 1090 in Figure 6b) were only variable at a single location and differed in allele frequencies across sites (i.e., depths) at that one location (Figure 6b). In contrast, other SNPs reflected potential interactions between depth and SST. Several SNPs had alleles that were fixed at either end of the latitudinal gradient but were variable in the middle and differed across depths at one of the mid‐range locations. In most cases, allele frequencies at the deep sites were similar to those further south while those in the shallow sites were similar to those from further north. However, this was not universal with clear exceptions that showed the opposite pattern (e.g., locus 6060 in Figure 6b). We did not identify any SNPs that had universally different allele frequencies across depths—for example, alleles that were consistently more frequent in deep environments across all locations.

Loci potentially under selection across environments. (a) Examples of loci significantly associated with maximum SST and reflecting underlying latitudinal structure. (b) Examples of depth‐associated SNPs.

Discussion

4

Our study provides numerous insights into the genetic diversity, connectivity and local adaptation of Eisenia arborea populations across both latitude and local‐scale gradients in depth. The results reveal clear population structuring driven by a combination of neutral and adaptive processes, shedding light on the microevolutionary responses of this foundation species to environmental variation at different spatial scales. Genetic diversity and demographic inferences drawn from the genomic data are partially consistent with predictions for marginal populations from the CMH (Eckert et al. 2008), with the lowest heterozygosity and highest inbreeding observed in the warm range‐edge populations. However, counter to expectation, the southern range edge populations had less differentiation than mid‐range populations that were of similar distance. We also demonstrated that despite high gene flow across sites (i.e., deep and shallow) in the same location, some SNPs did have different allele frequencies across depths that reflect putative loci under selection. However, these patterns were very location‐specific, suggesting that the role of depth in driving selection may vary across latitude or regional environmental contexts. Our results uncover patterns of population structure that have the potential to inform management of this species, for example through recognition of conservation units (CUs) based on genetic differentiation (Fraser and Bernatchez 2001; Coates et al. 2018) or in selection of appropriate donor populations for restoration and mariculture activities (Coleman, Wood, et al. 2020; Wood et al. 2021; Zertuche‐González et al. 2022).

Genetic Diversity and Central‐Marginal Patterns

4.1

Patterns of genetic diversity and demographics at the warm range‐edge were largely consistent with CMH expectations for marginal populations. In particular, the consistently lower heterozygosity and higher inbreeding coefficients in the southernmost locations are consistent with predictions of reduced genetic diversity due to smaller, more isolated populations, and potentially stronger directional selection in marginal populations (Lesica and Allendorf 1995; Jiang et al. 2015; Kitamura et al. 2020). Consistent with this hypothesis, historical population sizes were also generally lowest at Bahía Magdalena, the true southern range edge of this species. However, similar historical population sizes were also inferred for Danver's Point in British Columbia, which is technically central range for this species. We hypothesise that because Barkley Sound (the region containing Danver's and Taylor locations) is the furthest south population in British Columbia (i.e., north of the large distributional gap), E. arborea populations in this region may be marginal and have experienced similar processes to a true range edge population in terms of historical demographic dynamics. Consistent with this interpretation, populations from San Quintín exhibited comparatively stable inferred demographic histories, aligning with their higher genetic diversity and intermediate position within the species' range.

Our findings of lower heterozygosity and higher inbreeding coefficients in the southernmost populations mirror patterns seen in some other habitat‐forming seaweeds, such as the closely related Ecklonia radiata, where reduced diversity and elevated inbreeding were detected at multiple equatorward range edges along different coastlines (Vranken et al. 2021; Minne et al. 2025). However, some other studies have found no clear pattern associated with the range edge of species (e.g., Gierke et al. 2023) or even increased genetic diversity towards the warm range edge (Edgeloe et al. 2025). For example, the lowest genetic diversity and highest inbreeding in Nereocystis luetkeana is in the central range, presumably due to dynamics associated with the large coastal Salish Sea limiting gene flow and population size (Bemmels et al. 2025) as well as demographic processes associated with past glaciation (Bemmels et al. 2026). Thus, collectively these patterns suggest that marginal populations may not necessarily be restricted to range edges but instead, may occur across the range of species (Brown 1984).

However, despite showing demographic signatures consistent with marginal populations, the southern range‐edge populations exhibited unexpectedly low genetic differentiation from one another compared to some mid‐range populations. This apparent decoupling between demographic marginality and genetic connectivity at the southern range edge highlights an important nuance in applying CMH to marine systems. While CMH predicts reduced genetic diversity and increased differentiation at range margins, these outcomes do not all need to co‐occur. In Eisenia arborea , the southernmost populations appear demographically marginal, as evidenced by low heterozygosity, elevated inbreeding and reduced historical population sizes at Magdalena, yet they remain genetically similar to one another. This pattern is consistent with smaller, environmentally stressed populations embedded within a relatively continuous and well‐connected stretch of habitat, potentially facilitated by regional currents (e.g., Low Latitude California Current; LLCC). Under such conditions, episodic recruitment failures and strong selection at the warm edge may erode within‐population diversity without necessarily leading to spatial isolation from populations further north. These results suggest that demographic marginality and genetic isolation can be decoupled in marine foundation species, particularly where habitat continuity and alongshore transport persist at range edges.

While it has previously been hypothesised that British Columbia E. arborea populations may represent a recent introduction associated with shipping (Druehl 1970), our results do not show evidence of this. Metrics of genetic diversity were similar between British Columbia and Mexican samples, whereas we would have expected to see reduced genetic diversity in a recent non‐native introduction (Sax et al. 2007). Although our demographic analysis in StairwayPlot2 suggests that historical population declines have occurred across most sites over the recent generations (Figure 4), there was not a unique and more recent bottleneck at the Danver's Point shallow site (the only British Columbia site included in that analysis due to sample size limitations—see Section 2). Moreover, these bottlenecks are on the timescale of hundreds to thousands of years (albeit with great uncertainty) rather than on a timescale that would be expected to be driven by human‐induced dispersal. Instead, a more plausible hypothesis is that historical declines were driven by geological climate change, such as cycles of glaciation which would have made large parts of the northern range of this species uninhabitable and led to rapid cooling in southern parts of its range (Grant and Bringloe 2020; Mann and Gaglioti 2024). British Columbia populations were also the most genetically distinct of those sampled, making it clear that these populations are not recent invasions from the Baja California peninsula. We note that further sampling in California would help to better understand migration pathways across the range of E. arborea .

Local Adaptation Across Spatial Scales

4.2

As kelp forests continue to face threats from climate change and other anthropogenic stressors (Filbee‐Dexter and Wernberg 2018; Wernberg et al. 2019, 2024), understanding the genetic basis of adaptation across environmental gradients becomes increasingly important (Coleman, Wood, et al. 2020; Wood et al. 2021). In our study, the high connectivity across depths within locations indicates substantial vertical gene flow. This high gene flow raises questions as to whether depth‐related selective pressures might be insufficient to overcome the homogenising effects of gene flow at the local scale. However, our environmental association analyses (EAA) identified several depth‐associated, potentially adaptive loci, suggesting that selection is indeed affecting genetic structure across depths despite this high gene flow. For example, the RDA axis constrained to depth explained 1.3% of the variation (Figure S7) which, although small, is still consistent with some local adaptation across depths. Moreover, many of the loci associated with depth in our analyses exhibited highly site‐specific patterns of allele‐frequency divergence rather than consistent shifts across all locations. Perhaps this is because depth is a composite gradient encompassing multiple parameters (e.g., light, temperature, nutrient availability), suggesting that fine‐scale adaptation to depth may reflect interactions between local environmental conditions and broader climatic context, rather than a simple, generalisable depth effect.

We found much clearer patterns of differentiation with temperature, when compared to depth. Temperature‐associated loci generally showed consistent latitudinal transitions, albeit with variation in the specific allele transition points, suggesting that thermal adaptation is a key driver of genetic differentiation in E. arborea . This finding is especially relevant given that E. arborea is distributed further south than any other kelp species in North America (Setchell and Gardner 1925; Druehl 1970), implying particularly high thermal tolerance that may be underpinned by these temperature‐associated loci or genes linked to them (Vranken et al. 2021; Minne et al. 2025). The detection of these potentially adaptive loci is critical in the context of intensifying marine heatwaves and ocean warming, which are major threats to kelp forests worldwide (Wernberg et al. 2024), including in the Northeast Pacific (Rogers‐Bennett and Catton 2019; Beas‐Luna et al. 2020; Starko, Timmer, et al. 2024). While there was relatively little overlap between loci identified by LFMM and RDA, we chose to include loci detected by both methodologies based on a growing consensus in the field that EAA methods can detect different, yet biologically meaningful, signals of selection (Capblancq et al. 2018; Forester et al. 2018; Capblancq and Forester 2021). Disregarding loci identified by only one method risks omitting true positives and underestimating the complexity of local adaptation. Instead, we interpret the results of each method as providing complementary perspectives on selection. Nonetheless, we acknowledge that this approach may increase the false positive rate and encourage future studies to validate candidate loci functionally or experimentally.

Implications for Conservation and Management

4.3

Our findings have important implications for conservation and management of this important kelp species in the face of climate change. The observed genetic structure and genomic evidence for local adaptation suggest that E. arborea populations likely have different characteristics (e.g., thermal tolerance) and may therefore respond differently to environmental change (Coleman, Wood, et al. 2020). The limited connectivity inferred across locally adapted populations indicates that management strategies should consider distinct genetic units (i.e., CUs) particularly when planning restoration efforts (Coleman, Wood, et al. 2020; Wood et al. 2021) and may consider assisted gene flow approaches in the future if populations fail to track the rapid pace of global ocean warming (Coleman, Wood, et al. 2020). For example, translocating warm‐adapted genotypes from the southern part of the range of E. arborea to more northern areas may help bolster the resilience of mid‐range populations to climate change (by directly increasing their thermal tolerance) while also helping to preserve the uniquely adapted genotypes found towards the warm range‐edge of the species (Coleman, Wood, et al. 2020). However, it would be critical to experimentally validate these inferred functional differences between populations, for example using thermal tolerance assays, prior to implementing assisted gene flow strategies.

While southern populations may be locally adapted to warm conditions through strong selection and isolation, the combination of high homozygosity and restricted connectivity to populations further north suggests that populations may lack the genetic diversity needed to adapt further as temperatures rise beyond current physiological thresholds. This is especially concerning given projections of increased thermal stress in this region (Pörtner et al. 2019) and the low potential for natural gene flow from elsewhere in the range (see Figure 2). Although some small kelp populations may persist despite high inbreeding (Bemmels et al. 2025), elevated fixation likely reduces the scope for evolutionary responses (Lande and Shannon 1996; Willi et al. 2006). Overall, the convergence of historical bottlenecks, high inbreeding and limited connectivity emphasises the need for conservation strategies that account for both evolutionary history and spatial genetic structure when managing kelp forests under climate change.

Caveats and Future Directions

4.4

While our study provides valuable insights into the genetic structure and adaptive potential of E. arborea , several limitations should be acknowledged. The relatively small sample sizes at some locations, particularly for deep samples, may have limited our ability to detect fine‐scale patterns of genetic variation. Therefore, the fact that we successfully identified several depth‐associated loci despite these limited sample sizes speaks to the potential strength of selection to drive differentiation across even small spatial distances. Another limitation is that our EAA analyses focused only on maximum SST and depth, but other environmental factors such as light availability, nutrient concentrations and wave exposure are known to influence kelp adaptation (e.g., Coyer 2004; Starko et al. 2020). Moreover, we only looked at two different depths from each location and did not always sample the deepest populations, which can extend to 30 m depth along the Baja California peninsula (but in Barkley Sound only extends to about 10 m). This may have limited our ability to capture signals of local adaptation to depth.

Future research should expand sampling efforts to include a wider range of conditions both within and across locations, including a wider range of depths and other environmental variables. Increased sampling across the broad geographic range of this species, including California and northern British Columbia, would be highly informative and would provide insights into the historical and contemporary processes driving the disjunct distribution. For example, E. arborea is distributed throughout Southern California and in Monterey Bay but the relationship between populations in these areas and those included in our study remain unclear. Future research that integrates fine‐scale environmental data and physiological performance across both depth and latitudinal gradients will be particularly valuable for linking genotype to phenotype in this species. Genetic divergences should be examined alongside phenotypic plasticity analyses at various biological levels (e.g., metabolome, ecophysiology and morphology). Finally, specific experimental approaches would help to determine how different genotypes respond to environmental variables that influence their distribution, including ocean warming associated with climate change.

A key biological consideration is the biphasic life cycle of kelps, in which microscopic gametophytes may persist through periods that are unfavourable for sporophytes and can exhibit higher tolerance to thermal and other stressors than adult sporophytes (e.g., Muth et al. 2019). This ‘life‐stage buffering’ could reduce the extent to which standing sporophyte abundance reflects extinction risk and may complicate simple links between contemporary genetic diversity and near‐term resilience. However, because sporophytes represent the dominant habitat‐forming stage and because selection may act on traits expressed in either life stage, genomic structure and genotype‐environment associations observed in sporophytes remain informative for identifying differentiated units and potential signatures of local adaptation. We therefore interpret our results as describing adaptive and neutral structure in the standing sporophyte populations, while recognising that gametophyte dynamics may modulate demographic sensitivity to warming and merit further targeted study.

Conclusions

5

Our study reveals that both neutral and adaptive processes have shaped the genetic structure of Eisenia arborea populations across multiple spatial scales. The complex interplay between broad‐scale latitudinal gradients and fine‐scale depth‐related environmental variation has created a mosaic of genetic diversity and local adaptation that will likely influence how this foundation species responds to ongoing environmental change. The combination of limited connectivity at broad spatial scales, evidence for local adaptation to both latitudinal temperature gradients and depth, and consistently high inbreeding coefficients indicates that a complex interplay between gene flow, selection and genetic drift drives evolutionary processes in this species. Future work integrating functional assays or experimental validation would be valuable to corroborate the putatively adaptive loci identified here. Nonetheless, these findings provide crucial insights into the genetic architecture and adaptive potential of E. arborea , with important implications for predicting and managing the responses of these vital kelp forests to ocean warming. As marine environments continue to change rapidly, understanding the genetic basis of adaptation to environmental variation will be essential for developing effective conservation strategies to maintain marine ecosystem health into the Anthropocene.

Author Contributions

Samuel Starko: conceptualization (equal), data curation (equal), formal analysis (lead), investigation (equal), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Thomas Wernberg: conceptualization (equal), funding acquisition (lead), investigation (equal), writing – review and editing (supporting). Jose Miguel Sandoval Gil: conceptualization (equal), investigation (equal), project administration (equal), resources (equal), writing – review and editing (supporting). Jose Zertuche‐González: conceptualization (equal), investigation (equal), methodology (equal), project administration (equal), writing – review and editing (supporting). Ricardo Cruz‐López: conceptualization (equal), investigation (equal), methodology (equal), project administration (equal), writing – review and editing (supporting). David Wheeler: formal analysis (supporting), methodology (equal), writing – review and editing (supporting). Jacqueline Batley: investigation (supporting), project administration (supporting), resources (supporting), writing – review and editing (supporting). Melinda A. Coleman: conceptualization (equal), funding acquisition (equal), investigation (equal), methodology (equal), project administration (equal), resources (equal), supervision (equal), writing – review and editing (equal).

Funding

This work was supported by the Forrest Research Foundation, CONACYT‐Ciencia Básica Project (CLIMAVEMAR A1‐S‐8382) and Australian Research Council (DP200100201, DP240100230).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Table S1: Mean diversity indices associated with each sample. Figure S1: Histogram of error rates from duplicate samples (n = 8 sets). Most samples (n = 6; 75%) had an error rate below 1%. Figure S2: sNMF structure analysis with K = 5. Each bar indicates a different individual coloured by its inferred ancestral population. Figure S3: Cross entropy from sNMF structure analysis. Note the ‘elbow’ at K = 4. Figure S4: Isolation by distance plot showing a significant relationship between geographic distance and genetic differentiation. Each datapoint represents a pairwise comparison of two sampling sites. Figure S5: Reconstruction of effective population size (N e) through evolutionary time using EPOS. Note that while the MBS population (a) appears to have a larger population size at present, the model error includes a value as low as 1 for the most recent generation, which would represent an instantaneous decline from the lower end of the confidence bands at timepoint 1. Note that the x‐axis indicates the number of generations into the past. Figure S6: Venn diagram showing SNPs detected as putatively under selection. (a) SNPs associated with environmental variables broken down by driver and method. (b) Venn diagram showing SNPs identified using all EAA methods and BayeScan for outlier detection. Figure S7: Redundancy analysis (RDA) with axes constrained to sea surface temperature and depth. Each point represents an individual where colour and shape indicate the sampling site and depth, respectively.