Local thermal adaptation and limited gene flow constrain future climate responses of a marine ecosystem engineer

Abstract Rising ocean temperatures and extreme temperature events have precipitated declines and local extinctions in many marine species globally, but patterns of loss are often uneven across species ranges for reasons that are poorly understood. Knowledge of the extent of local adaptation and gene flow may explain such patterns and help predict future trajectories under scenarios of climate change. We test the extent to which local differentiation in thermal tolerance is influenced by gene flow and local adaptation using a widely distributed intertidal seaweed (Hormosira banksii) from temperate Australia. Population surveys across ~2,000 km of the species range revealed strong genetic structuring at regional and local scales (global F ST = 0.243) reflecting extremely limited gene flow, while common garden experiments (14‐day exposures to 15, 18, 21°C) revealed strong site differences in early development and mortality in response to elevated temperature. Embryos from many sites spanning a longitudinal thermal gradient showed suppressed development and increased mortality to elevated water temperatures, but populations originating from warmer and more variable thermal environments tended to be less susceptible to warming. Notably, there was significant local‐scale variation in the thermal responses of embryos within regions which was corroborated by the finding of small‐scale genetic differences. We expect the observed genetic and phenotypic differentiation to lead to uneven responses to warming sea surface temperatures in this important marine foundation species. The study highlights the challenges of predicting species responses to thermal stress and the importance of management strategies that incorporate evolutionary potential for “climate‐proofing” marine ecosystems.


| INTRODUC TI ON
The speed and magnitude of biotic shifts being triggered by climate change pose a major challenge for marine biodiversity conservation.
Rising ocean temperatures, increasing acidification, and changing ocean currents are contributing to fundamental and irreversible ecological transformations in marine ecosystems at a global scale (Babcock et al., 2019;Harris et al., 2018;Hoegh-Guldberg & Bruno, 2010). Species living close to their physiological limits are of particular concern and will become increasingly dependent on their ability to overcome environmental change via dispersal, physical tolerance, and evolutionary adaptation Portner & Gutt, 2016). Understanding the environmental resilience of ecosystem engineers (i.e., keystone and foundation species) is of particular importance, given their response to environmental change will have major impacts on community structure and ecosystem function (Babcock et al., 2019;Hoegh-Guldberg & Bruno, 2010; 2019; Thomson et al., 2015). Such information is necessary to inform marine conservation planning directed at minimizing biodiversity loss, and preserving environmental, economic, and cultural values (Magris, Pressey, Weeks, & Ban, 2014).
Macrophytes are key foundation species in temperate oceans underpinning reef ecosystem services (Dayton, 1972;Jones, Lawton, & Shachak, 1994) but are at increasing threat from environmental change with declines seen as a result of warming (Vergés et al., 2016(Vergés et al., , 2014, heatwaves (Thomson et al., 2015;Wernberg, Bennett, et al., 2016), over-harvesting (Steneck et al., 2002), and urbanization (Airoldi & Beck, 2007;Coleman, Kelaher, Steinberg, & Millar, 2008;Connell et al., 2008). Critically, such declines appear to be long term or even permanent, highlighting a need to understand the adaptability and vulnerability of marine macrophytes to future change. This is pertinent in areas where marine macrophyte communities are showing signs of acute and chronic climate stress, such as climate change "hot spots" that are particularly prone to extreme temperature events, and that are experiencing sea surface temperature orders of magnitude above the global average (Babcock et al., 2019;Hobday & Pecl, 2014;Wernberg et al., 2011). This is a major concern as thermal stress is known to impair photosynthetic, respiratory, and cellular function in macrophyte species, inhibiting growth, inducing mortality and disease, and resulting in range contractions and local extirpations (Flukes, Wright, & Johnson, 2015;Mathieson & Dawes, 1986;Pineiro-Corbeira, Barreiro, Cremades, & Arenas, 2018;Smale & Wernberg, 2013;Wernberg, Bettignies, Joy, & Finnegan, 2016). Furthermore, thermal stress can suppress macrophyte resilience to environmental stressors associated with other anthropogenic and natural perturbations (Wernberg et al., 2010) as well as alter key ecological interactions that impact macrophyte populations, such as grazing and predation (Miranda et al., 2019;Provost et al., 2017;Vergés et al., 2014). Coupled climateecosystem models suggest that declines are likely to intensify in coming decades, with projections of major range contractions of temperate seaweeds, and potential risks of extinction for many species (Martinez et al., 2018).
The resilience of marine macrophyte populations to environmental change will be largely determined by levels of standing adaptive genetic variation and patterns of gene flow among populations.
Traditionally, insights into climate resilience have been gained by characterizing patterns of gene flow among populations distributed across thermal gradients. While gene flow can be an impediment to local adaptation, it can also assist the adaptation process by making thermally adapted genotypes available for selection Hoffmann & Willi, 2008). Also, quantitative approaches (e.g., common garden and transplant experiments) are used to assess the extent to which species and populations have adapted to different environments, by testing the genetic basis of trait variation spanning environmental gradients (Hoffmann & Willi, 2008).
Such studies have been particularly useful for assessing historical responses to environmental change and predicting the evolvability of species from both marine and terrestrial systems (Sgro, Lowe, & Hoffmann, 2011;Sherman & Ayre, 2008), but they are limited for marine macrophytes.
Predicting species responses to warming environments requires an understanding of intraspecific variation in thermal responses and gene flow across species ranges (King, McKeown, Smale, & Moore, 2018). Predictive tools such as species distribution and climate niche models do not directly consider physiological variation unless this variation contributes to extant distributions (Razgour et al., 2019;Willis et al., 2015). Such models may therefore be limited in predicting range shifts, particularly at edges of a species' distribution, where temperatures are expected to exceed theoretical species "thermal niche" limits or safety margins (Bennett, Wernberg, Joy, Bettignies, & Campbell, 2015), which may be countered by evolutionary shifts (Bush et al., 2016).
Variability in thermal physiology between populations has been established for several marine macrophytes (reviewed in King et al., 2018), which may contribute to uneven population responses to thermal stress (Bennett et al., 2015;Carballo, Olabarria, & Osuna, 2002;Filbee-Dexter, Feehan, & Scheibling, 2016;Starko et al., 2019;Tegner & Dayton, 1987;Thomsen et al., 2019). This variation often coincides with fine-scale genetic structuring (King et al., 2018), but empirical tests of this are lacking. Information on the genetic structure of populations and phenotypic differences among populations is needed to help understand the uneven thermal response of populations, identify those most vulnerable to thermal stress, and indicate where genotypes resilient to future climate changes might reside. These are all important components of climate change adaptation and restoration programs (Foden et al., 2019;Jordan, Hoffmann, Dillon, & Prober, 2017;Willis et al., 2015;Wood et al., 2019).
Hormosira banksii (Turner) Decaisne is the dominant intertidal macrophyte across temperate Australasia and is an autogenic ecosystem engineer with no functional equivalents (Bishop et al., 2009;Pocklington, Keough, O'Hara, & Bellgrove, 2019;Schiel, 2006). This species is highly sensitive to environmental disturbance associated with coastal development and urbanization (Doblin & Clayton, 1995;Keough & Quinn, 1998;Kevekordes, 2000), and thermal stress (Alestra & Schiel, 2015;Tait & Schiel, 2013). H. banksii embryo development is particularly sensitive to upwards shifts in temperature, with water 3°C higher than ambient levels causing significant mortality (Alestra & Schiel, 2015;Clark, Poore, Ralph, & Doblin, 2013). Yet significant genotype × environment interactions for embryo growth have been reported, as well as growth and photosystem traits that show heritable genetic variation relating to temperature (Clark, Poore, & Doblin, 2018;Clark et al., 2013). These findings suggest there is potential for selection to result in the evolution of more thermally tolerant populations. However, despite potential standing genetic variation in some populations for thermal adaptation, gene flow may be limited in H. banksii, restricting the adaptive potential of this species. Previous studies indicate strong genetic structuring among populations from the central region of Australia's east coast (Coleman et al., 2011;Coleman, Clark, Doblin, Bishop, & Kelaher, 2019) suggesting that gene flow is unlikely to facilitate the natural movement of genetic variants across temperature gradients to aid adaptation and recovery of depleted populations. Such genetic structuring could also contribute to local adaptation and variation in thermal tolerances among populations spanning thermal gradients in Australian temperate waters. However, the genetic structure across the full range of this species, and levels of population variability in thermal tolerance, is currently unknown.
In this study, we assessed the potential adaptability of H. banksii to rising ocean temperatures across much of its distribution. First, we use microsatellite markers to investigate patterns of genetic diversity, gene flow, and population connectivity from sites spanning ~2,000 km of the distribution of this species around mainland Australia and covering a wide temperature gradient (mean summer sea surface temperatures ranging from 16 to 24°C). Second, we undertake a common garden experiment to explore spatial variation in thermal responses among populations with different thermal histories by comparing embryo development at different experimental temperatures. Outcomes from these analyses will help determine the role of gene flow and local adaptation in affecting responses of populations to thermal stress. We expected populations under thermal extremes in the sampling range to be phenotypically adapted to those temperature environments, particularly in the presence of low levels of gene flow. Any adaptive changes would produce uneven thermal responses among populations across the species range. We discuss the implications of our results in the context of the future evolutionary potential of H. banksii, and conservation/ restoration approaches that might be considered to catalyze adaptation and recovery processes in light of climate change.

| Study species
Hormosira banksii is a perennial, dioecious, habitat-forming fucoid macroalga with a distribution extending >2,000 km from Albany in Western Australia to Lennox Head in New South Wales on mainland Australia, around Tasmania, the North and South Islands of New Zealand, and some of the smaller offshore islands in southern Australasia (Clark et al., 2018;Osborn, 1948). H. banksii is fertile throughout the year, releasing gametes on emersion during low tide (Levring, 1949). Once fertilized, eggs are thought to attach rapidly, thereby limiting dispersal (Bellgrove, Clayton, & Quinn, 1997. The thallus of H. banksii consists of multiple chains of hollow, water-filled, vesicles that arise from a diffuse holdfast (Osborn, 1948) capable of vegetative regeneration (Keough & Quinn, 1998;Schiel & Taylor, 1999), with individual plants surviving >11 years (Kain, 2015). H. banksii shows significant morphological variability among populations and across geographic regions, and between marine and estuarine environments shown to have functional consequences (Clark et al., 2018). This is thought to be correlated with local envi-

| Microsatellite analysis
A total of 373 H. banksii individuals from 26 sites spanning New South Wales, Victoria, and South Australia were genotyped at 10 microsatellite loci previously developed by Bellgrove et al. (2017). Loci were amplified by multiplex PCR using Eppendorf Mastercycler S gradient units following the protocol described by Blacket, Robin, Good, Lee, and Miller (2012). Genotyping was performed using an Applied Biosystems 3730 capillary analyzer (AGRF, Melbourne Australia), and microsatellite profiles were examined and scored manually using GENEMAPPER version 4.0 (Applied Biosystems). Independence of loci (absence of linkage) was assessed using the software GENEPOP, version 3.4 (Raymond & Rousset, 1995) and the probability test function, which found no significant association among loci. The software MICRO-CHECKER (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004) was used to assess microsatellite loci for null alleles and scoring errors using formula 1 outlined by Brookfield (1996) and found no evidence of null alleles across sites and loci. Descriptive statistics were calculated for the microsatellite data using FSTAT version 2.9.3 (Goudet, 1995), including allelic richness per population averaged over loci, Weir and Cockerham's inbreeding coefficient (F IS : the deficiency of heterozygotes relative to the level expected with random mating), and a global estimate of population differentiation (F ST ) with 95% confidence limits (Weir & Cockerham, 1984). In addition, population pairwise measures of F ST and their significance were determined using permutation (10,000 permutations). In order to overcome potential limitations of F ST calculations using multiallelic loci (Jost, 2008), additional estimates of population differentiation, global D est , and population pairwise measures of D est (significance determined using 10,000 permutations) were generated using GenAlEx 6.5 (Peakall & Smouse, 2006). The false discovery rate (FDR) procedure (Benjamini & Hochberg, 1995) was used to adjust significance levels when performing multiple simultaneous comparisons.
Estimates of observed (Ho) and expected (He) heterozygosity were determined using the Excel Microsatellite Toolkit (Park, 2001) and deviations from Hardy-Weinberg equilibrium (HWE) were determined using Genepop version 3.4 (Raymond & Rousset, 1995). An analysis of molecular variation (AMOVA) was performed in GenAlEx 6.5 (Peakall & Smouse, 2006) using pairwise F ST as the distance measure, with 10 000 permutations and missing data for loci set at 10%. Bayesian analyses implemented in STRUCTURE (Pritchard, Stephens, & Donnelly, 2000) were conducted to estimate the number of populations within the sample data. STRUCTURE identifies the number of distinct population clusters, assigns individuals to clusters, and identifies migrants and admixed individuals using genetic data only. To determine the number of populations (K), five independent simulations for K = 1-18 with 100,000 burn-in and 1,000,000 data iterations were run. Analyses were performed using the admixture model of population structure (i.e., each individual draws some fraction of their genome from each of K populations) with no prior population information included, and allele frequencies set as independent among populations. The most likely K was estimated using Evanno's ΔK (Evanno, Regnaut, & Goudet, 2005) in Structure Harvester (Earl & Vonholdt, 2012).
Rates of recent migration were estimated among each of the sampling locations using a Bayesian algorithm implemented in BAYESASS v.3.0.3 (Wilson & Rannala, 2003). The program estimates migration among populations within the last three generations. To identify movements among populations, five independent runs of 10 7 Markov chain Monte Carlo (MCMC) iterations were used following a burn-in period of 10 7 and a sampling interval of 500 steps.
Chains were compared with a stationary posterior distribution for  Figure 1c). Sea surface temperature rather than air temperature was F I G U R E 1 Maps showing Hormosira banksii sampling locations for the population genetic analysis (map "B") and common garden experiments (map "C," modified from Ierodiaconou et al., 2018). Map c shows the nesting of collection sites within three distinctive temperature zones (low, intermediate, and high), with the coastline color coded according to local mean summer sea surface temperatures (SST). Summary plots of the population subdivision based on STRUCTURE analyses of Hormosira banksii multilocus genotypes are also included ("d" and "e"). Single vertical lines broken into segments reflect each individual in the STRUCTURE summary plot, where segments are proportional to the membership coefficient for each of the population clusters. Individuals are arranged into sites from which they were sampled following the order given in Table 1 (site codes are also derived from Table 1). STRUCTURE plots are provided for the first hierarchical level of structure (all states included in the analysis; plot "d") and the second hierarchical level of structure (independent analysis of sites from New South Wales, Victoria, and South Australia; plot "e") Phillip Bay were used only for the common garden experiment.

| Thermal response experiments
H. banksii zygotes from each of the three defined thermo-geographic regions (Southwest, East, and PPB) were reared under common culture conditions for experimental purposes following the collections outlined in Section 2.3.1. Thermal response experiments were performed by exposing zygotes from each site to three controlled temperature treatments, 15, 18, and 21℃. These temperatures were based on 15 and 18℃ falling within the bounds of annual sea surface temperatures (SST) across all regions, while only sites from the PPB region are exposed to average summer SSTs as high as 21℃ (Table 2).
Initially, gamete release for all sample sites was stimulated in accordance with Kevekordes and Clayton (1996), whereby thalli were rinsed in freshwater, patted dry using paper towel, and exposed to ambient temperature and light for up to 1 hr until gametes were exuded.
Eggs were collected prior to sperm to maximize fertilization success, as sperm viability is known to deteriorate quickly upon release into seawater (Doblin & Clayton, 1995;Levring, 1949 (Kevekordes & Clayton, 1996). Once standardized, gamete solutions from each source population were mixed to initiate fertilization (within population crosses only), with 700 µl from each zygote solution immediately pipetted onto 6 individual glass coverslips placed inside each of 18 replicate Petri dishes (6 for each temperature treatment) and left to settle for approximately 5 min (Kevekordes & Clayton, 1996). Twenty mL of 15℃ microfiltered seawater was gently added to each Petri dish to submerge the glass coverslips, after which six Petri dishes from each sample location were randomly assigned to one of three tempera- TA B L E 2 Location details for 9 collection sites of Hormosira banksii from southeastern Australia used for the common garden experiment, including mean summer and mean seasonal ranges in sea surface temperature in decimal degrees (°C) insufficient zygotes had settled onto a single coverslip for scoring, a second one was sampled. At day 7, photographs of embryos were taken with a mounted camera (ZEISS-Axiocam ERc 5s) at ×100 magnification and subsequently 5 embryos randomly selected from each replicate for measurement of the rhizoid lengths using the ImageJ analysis package (version: 2.0.0-rc-43/1.51p, https ://imagej.net).

| Statistical analyses for common garden experiment
Fertilization success prior to temperature manipulation was analyzed by a two-factor nested analysis of variance (ANOVA) with re- , where σ 2 b and σ 2 w are the respective phenotypic variances between and within groups of populations, c is an estimate of the proportion of the total variance due to additive genetic effects across populations, and h 2 is narrow-sense heritability, the proportion of phenotypic variance due to additive genetic effects (Brommer, 2011). The c/h 2 default level of 1 was used; 95% confidence intervals were calculated with 1,000 bootstrapping data frames. As values of c/h 2 are often not known for wild populations and can be variable among populations, the robustness of this ratio on phenotypic divergence was evaluated using additional c/h 2 values of 0.5 and 0.1. As microsatellite genotypes were not available for sites included in the common garden experiment, we contrasted P ST against global F ST , and against estimates of genetic divergence between and within regions.

| Population genetic analysis
Across the 10 microsatellite loci and 18 sample locations, we detected a total of fifty-one alleles (one of which was exclusive to South Australian sites), with a mean of 6.0 alleles per locus (range 4-8). Levels of genetic diversity (allelic richness and heterozygosity) were largely consistent across sites over all loci and heterozygosity, with a mean allelic richness of 2.59 (range 2.16-3.00) and a mean observed heterozygosity of 0.41 (range 0.31-0.48; Table S1). Most sites were found to conform to HWE (Table S1) (Table 3).
Similarly, AMOVA analyses indicated that a significant proportion of the genetic variance (9%) could be attributed to difference among regions (p < .01, F ST = 0.092), and 16% of the variance due to differences among sites within regions (p < .01, F ST = 0.166). The relationship between genetic and geographic distance ( Figure S1 Genetic differentiation between regions and sites within regions is depicted by the discriminant analysis of principal components (DAPC) of the microsatellite variation ( Figure S2). When all sites were included in the analysis, clear regional differences were evident, with sites from NSW, VIC, and SA separating across the axes and representing three distinct clusters (estimated using the find cluster function; Figure S2a). DAPCs performed on sites from each state show regional genetic structuring with sites from common regions typically clustering more closely, but with centroids of population clusters rarely overlapping ( Figure S2a-d).  Table 1. STRUCTURE Bayesian clustering analyses confirmed significant genetic differentiation between the NSW, VIC, and SA sites, with analyses indicating assignment of individuals from the respective states to three distinct population clusters (ΔK = 3; Figure 1d). STRUCTURE analysis of each region suggested significant structuring with each of the regions, consistent with the DAPC analysis. Further structuring at the state scale was evident when STRUCTURE analyses were performed on sites from each region independently (Figure 1e). Analyses of NSW sample locations revealed three distinct population clusters (ΔK = 3) corresponding with the Central Coast, Sydney, and South Coast regions. Similarly, significant regional differences were observed in SA where individuals were assigned to three population clusters (ΔK = 3) corresponding with the Fleurieu Peninsula, Yorke Peninsula, and Eyre Peninsula regions. Regional genetic structuring was less evident in Victoria where two population clusters were determined (ΔK = 2), but not reflecting an obvious geographic pattern. We expect this is likely due to comparatively lower level of genetic differentiation among Victorian sample locations, and the limited ability of STRUCTURE to resolve fine-scale patterns of genetic structure (Evanno et al., 2005).

TA B L E 3 Pairwise estimates of
Results from the Bayesian migration analyses indicated limited migration among sites within the last three generations (Table S1).
Estimates of the strength and directionality of migration demonstrated that each site has been largely dependent on recruitment from local sources, with minimal migration from nonlocal sources.
Evidence of weak recent migration (average 8.3% migrants per generation) was recorded between three pairs of neighboring sites. In the Sydney region, evidence of unidirectional migration of 5% was found between LOR and SUT, while bidirectional migration of 9%-11% between sites WIL and ALD on the Fleurieu Peninsula separated by <2 km was also detected.

| Fertilization success
Fertilization success was high for all sites (mean ± SE: 92.7 ± 0.98%, N = 81 pooling across sites; Figure S3), and there was a marginally nonsignificant difference in fertilization success among regions (2-factor nested ANOVA: F (2,6) = 5.08; p = .051) and no significant difference among sites within regions (F (6,72) = 1.09; p = .377) at the beginning of the experiment. This was expected as all gamete mixes were initially prepared under common temperatures (15°C), where fertilization was expected to be largely complete within minutes (Kevekordes & Clayton, 1996), prior to exposure to experimental temperatures.

| Embryo development and rhizoid growth (Day 7)
The response of developing embryos to different temperature treatments was complex and variable among sites and regions. The percentage of mature embryos at day 7 (>4 cell divisions and presence of rhizoids) among temperature treatments were consistent across sites within regions (ANOVA temperature × site(region); F (12,135) = 0.77; p = .685) but not among regions (ANOVA temperature × region; F (4,12) = 8.58; p = .002). There was no significant difference in the percentage of mature embryos among temperature treatments for sites from the PPB region (SME temperature F (2,12) = 2.352; p = .137; Figure 2a), although this region did have the lowest proportion of mature embryos compared with sites from the Southwest and East regions. In contrast, the percentage of mature embryos differed between temperatures for the F I G U R E 2 Development of embryos of Hormosira banksii from three source regions (southwest Victoria, eastern Victoria/southern NSW, and Port Phillip Bay, pooling three sites in each region) after 7 days of culture at 15, 18, and 21°C showing a) mean (±SE) percentage of embryos with greater than four cell divisions and well-developed rhizoids, and (b) mean (±SE) rhizoid length. Broken horizontal lines above bars in (b) indicate significant differences among temperature treatments based on Tukey's test after the full ANOVA model, whereas letters on bars in (a) indicate Tukey's test results among treatments for each region individually after simple main effects analysis (because of significant temperature × region interaction in the full model) Southwest and East regions, with significantly lower development in the warmest treatment compared with 15 and 18°C (Tukey's HSD, p < .05; Figure 2a). Rhizoid lengths increased significantly with temperature, with length being greatest at 21°C (ANOVA temperature; F (2,4) = 13.91; p = .001; Tukey's HSD p < .05 for all pairwise comparisons); this pattern was consistent across sites within regions (ANOVA temperature × site(region); F (12,134) = 1.66; p = .082) and among regions (ANOVA temperature × region interaction; F (4,12) = 2.35; p = .114; Figure 2b).

| Embryo mortality and meristematic activity (Day 14)
Embryo mortality by day 14 averaged 12.3 ± 0.84% (N = 162) but varied among sites within regions with temperature (ANOVA temperature × site(region) interaction; F (12,135) = 3.12; p = .001; Figure 3), which was largely driven by significantly higher mortality at 21°C than either 15 or 18°C for one site (Backyards) in the Southwest (Tukey's HSD, Bonferroni-p = .017 & .004, respectively; Figure 3). Of the embryos that survived to day 14, meristematic activity (the presence of apical hairs) among temperature treatments was not consistent among sites within regions (ANOVA temperature × site(region); F (12,135) = 2.66; p = .003; Figure 5), but temperature effects on meristematic activity were consistent among regions (ANOVA temperature × region; F (4,12) = 2.55; p = .093). Overall, there were significantly more embryos with apical hairs at day 14 in the 21°C treatment compared with both 18 and 15°C treatments (Tukey's HSD, p < .001 for 21 versus 18°C and 21 versus 15°C), but this pattern was only significant for 2 sites in each of the PPB and East regions (Tukey's HSD after SME temperature; Figure 4).

| Regression analyses and P ST -F ST comparisons
Regression analyses suggest an association between home site temperature environments and thermal response in the common garden experiment. We found a significant negative relationship between home site mean seasonal variation in sea surface temperature and 7-day rhizoid development (r = −.64, p = .007) and 14-day mortality (r = −.65, p < .001) at 21°C ( Figure 5). Nonsignificant negative trends were also noted when home site mean summer SST was regressed against 7-day rhizoid development (r = −.45, p = .225) and 14-day F I G U R E 3 Mean (±SE) percentage mortality of embryos of Hormosira banksii from three random sites in each of three source regions (southwest Victoria, eastern Victoria/southern NSW, and Port Phillip Bay) after 14 days of culture at 15, 18, and 21°C. Letters above bars indicate Tukey's test results among temperature treatments for each site individually after simple main effects analysis (because of significant temperature × site(region) interaction in the full model) mortality (r = −.58, p = .104) at 21°C (Figure 5). There was no association between home site mean seasonal variation in sea surface temperature or site mean summer temperature with 7-day rhizoid length (r = .05, p = .905; r = .05, p = .908; respectively), and 14-day apical hair development (r = .18, p = .647; r = .10, p = .801; respectively) at 21°C.

| Population connectivity
Genetic and phenotypic differences among populations may contribute to uneven species responses to climate change; such unevenness may result from the genetic history of populations as well as selection associated with prior thermal history. Understanding the interplay of these factors helps predict thermal responses of species to increasing incidences of acute and chronic thermal stress. Within this context, our study on the dominant intertidal macrophyte in Australasia, H. banksii, indicates highly restricted gene flow across much of its contemporary distribution with evidence of strong genetic structuring on all spatial scales examined. These findings are consistent with those reported previously for H. banksii over smaller geographic ranges (Bellgrove et al., 2017;Coleman et al., 2011Coleman et al., , 2019 and reports of finescale genetic structuring in other marine macrophytes including fucoids (Coleman & Brawley, 2005;Williams & Difiori, 1996). Common garden experiments revealed a heterogeneous spatial distribution of thermally tolerant phenotypes, which may be related to local-scale thermal adaptation. These findings suggest that the impacts on this important foundation species of warming ocean temperatures under climate change may be uneven. Opportunities for recovery of depleted populations and adaptation in populations vulnerable to thermal stress through natural gene flow may be limited, but interventions such as translocation and assisted migration could assist the recovery or bolster the resilience of existing populations given the existence of relatively tolerant populations across the species range (Marzinelli, Leong, Campbell, Steinberg, & Verges, 2016;Wood et al., 2019).

| Heterogeneity in thermal response
Consistent with Clark et al. (2013) and Alestra et al. (2015), our study suggests the potential for selection to increase thermal tolerance across the species range. Our common garden experiments indicate strong differences in early development between sites sampled from three geographical regions in response to differences in temperature. High water temperatures (21°C) often led to suppressed embryonic development and increased mortality in individuals from sites from cooler regions, while others from warmer and more variable temperature regions were relatively unaffected. Our correlation analyses indicated negative linear relationships between homesite, sea surface temperatures (mean summer temperature and seasonal variation), and thermal response at 21°C. These findings suggest genetically determined variation in thermal tolerance across broad thermal gradients in H. banksii, consistent with results reported for other fucoid and laminarian macrophytes (Bennett et al., 2015;King et al., 2018;Mohring, Wernberg, Wright, Connell, & Russell, 2014;Staehr & Wernberg, 2009).
Our results also indicate site variation within regions in the thermal responses of embryos to temperature stress, reflecting potential adaptation on local spatial scales. Although there are few studies on micro-scale variation in marine macrophytes, signatures of salinity-associated adaptation among fucoid (Fucus serratus) populations separated by as little as 12 km have been previously demonstrated (Coyer et al., 2011), and reciprocal transplants of genotypes from low to high intertidal zones (10s of meters scale) in another fucoid (Silvetia compressa) demonstrated clear home site advantages (Hays, 2007). These findings add to a growing literature supporting the notion that selection processes operate on finer spatial scales than currently assumed in marine ecosystems, with local habitat heterogeneity contributing to genetic adaptation in a range of marine organisms (Babin, Gagnaire, Pavey, & Bernatchez, 2017;Miller et al., 2019;Sherman & Ayre, 2008).
While the majority of studies reporting phenotypic variation across environmental gradients in marine macrophytes have not considered the relative contributions of local adaptation and plasticity (King et al., 2018) instance, macroalgae often show plastic responses in morphology and physiology to changes in temperature (Flukes et al., 2015;Hurd et al., 2014;Reusch, 2014). However, in our case we provide evidence for adaptive heritable variation that can contribute to phenotypic variation alongside this plasticity. The development of heritable differences among sites may reflect limited gene flow in the species. In the absence of gene flow, populations are expected to evolve phenotypes in response to local selection pressures, helping to produce the patterns observed here. Local adaptation over short distances has been shown to occur in terrestrial plant systems when gene flow patterns are conducive to adaptation (e.g., Buehler et al., 2013;Byars, Papst, & Hoffmann, 2007).
Although we have interpreted differences among sites as reflecting genetic variation, such comparisons should ideally involve additional generations because of the possibility of cross-generation effects induced by environmental conditions (Schiffer, Hangartner, & Hoffmann, 2013). Such effects can lead to phenotypic changes spanning multiple generations (i.e., transgenerational plasticity), including those due to epigenetic mechanisms (Walsh et al., 2016), in the absence of fixed genetic differences that reflect changes in allele frequencies. Cross-generation studies will help determine the rela-  (Bennett et al., 2015;Carballo et al., 2002;Filbee-Dexter et al., 2016;Starko et al., 2019;Tegner & Dayton, 1987;Thomsen et al., 2019). Such variation is often explained by heterogeneous habitat features (i.e., exposure, depth, bathymetry, and reef geomorphology) that contribute to different thermal environments at local and regional scales (Ierodiaconou et al., 2018). However, habitat differences are also expected to contribute to genetic and phenotypic differences influencing the climate niche of local populations due to selection, particularly for species with low levels of dispersal and strong genetic structuring (King et al., 2018;Miller et al., 2019). This highlights the need for improved evolutionary information on marine macrophytes at the population level to help understand what drives differential responses among populations to thermal stress. Also, this highlights the importance of species distribution and climate niche models that account for population differentiation in order to generate more reliable predictions of species responses to future climates (Martinez et al., 2018;Razgour et al., 2018;Valladares et al., 2014).
Genomic sequencing approaches (i.e., reduced genome representation, whole genome, and transcriptome sequencing) will help in further elucidating patterns of differentiation across sites. Such methods provide unprecedented sensitivity for resolving fine-scale genetic structure, as well as identifying signatures of selection and putative candidate genes that underlie adaptive differences between species populations (Jordan et al., 2017;Miller et al., 2019).
Unfortunately, this study was restricted to the use of microsatellite markers due to DNA not being of sufficient quality for high-throughput sequencing methods. However, the 10 polymorphic loci used in this study sufficiently resolved patterns of population genetic structure, and microsatellite markers have proven effective for such analyses in other marine macrophyte systems (Reynolds et al., 2017;Wernberg et al., 2018).

| Assisted migration and adaptation
Our findings have implications for predicting the recovery potential of depleted H. banksii in areas affected by coastal development, urbanization, and thermal stress (Bellgrove et al., 2010;Keough & Quinn, 1998). The apparent lack of gene flow among sites within regions suggests opportunities for natural recolonization following local mortality are low (Coleman et al., 2008). Instead, interventions such as translocation and reseeding may be needed to assist recovery (Bellgrove et al., 2010;Campbell, Marzinelli, Verges, Coleman, & Steinberg, 2014). Moreover, strong population genetic structuring in H. banksii also suggests that gene flow is unlikely to assist populations in adapting to warming sea surface temperatures via the migration of thermally adapted genotypes (Alestra & Schiel, 2015;Clark et al., 2013). Instead, a lack of gene flow will necessitate in situ adaptation. While empirical studies suggest that maladaptation can be mitigated in the presence of weak gene flow across environmental gradients Sgro et al., 2011), our results suggest that little to no migration is occurring between populations separated by 10s of kilometers. Consequently, the future resilience of local populations to rising sea surface temperatures will depend on standing genetic variation and plasticity, in the absence of any intervention.
While our study suggests that thermally resistant genotypes can occur, selection may struggle to keep pace with a rapidly shifting climatic envelope and increasingly frequent extreme temperature events . This is particularly pertinent for macrophyte species in climate change hot spots such as southwestern and southeastern Australia. For locally adapted species with limited dispersal ability such as H. banksii, adaptive management strategies might be needed. These could include assisted migration of thermally adapted genotypes to populations showing signs of climate stress. Such approaches are being widely advocated as a tool for "climate proofing" threatened marine and terrestrial animal and plant communities (Aitken & Whitlock, 2013;Prober et al., 2015;Sgro et al., 2011). For species showing genetically determined clinal variation in thermal response, such approaches might involve the translocation of genotypes from warm environments into cooler areas (Bansal, Harrington, Gould, St, & B., 2015;Schueler et al., 2013).
However, alternative and more tailored approaches may be needed for species such as H. banksii where thermal adaptations across the species range are heterogeneous. In such cases, assisted migration strategies may require composite provenancing approaches, involving mixed genotypes from multiple source populations (Broadhurst et al., 2008;Prober et al., 2015;Weeks et al., 2011). This can help to broaden the genetic basis of introduced genotypes, providing at least some genetic variants that are preadapted to warmer oceanographic conditions. is acknowledged for early contributions to assessing the population genetics of H. banksii. GP Quinn is thanked for advice relating to statistical analyses, and M Young for assistance with geospatial aspects of the project.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Genetic and phenotypic data will be made publicly available in the DRYAD: https ://doi.org/10.5061/dryad.g4f4q rfmb (Miller, 2019).