Anthropogenic climate change has caused widespread loss of species biodiversity and ecosystem productivity across the globe, particularly on tropical coral reefs. Predicting the future vulnerability of reef-building corals, the foundation species of coral reef ecosystems, is crucial for cost-effective conservation planning in the Anthropocene. In this study, we combine regional population genetic connectivity and seascape analyses to explore patterns of genetic offset (the mismatch of gene-environmental associations under future climate conditions) in Acropora digitifera across 12 degrees of latitude in Western Australia. Our data revealed a pattern of restricted gene flow and limited genetic connectivity among geographically distant reef systems. Environmental association analyses identified a suite of loci strongly associated with the regional temperature variation. These loci helped forecasting future genetic offset in random forest and generalised dissimilarity models. These analyses predicted pronounced differences in the response of different reef systems in Western Australia to rising temperatures. Under the most optimistic future warming predictions (RCP 2.6), we observed a general pattern of increasing genetic offset with latitude. Under the most extreme climate scenario (RCP 8.5 in 2090-2100), coral populations at the Ningaloo World Heritage Area were predicted to experience a higher mismatch in genetic composition, compared to populations in the inshore Kimberley region. The study suggest complex and spatially heterogeneous patterns of climate-change vulnerability in coral populations across Western Australia, reinforcing the notion that regionally tailored conservation efforts will be most effective at managing coral reef resilience into the future.

1. DArT sequencing

Population samples were collected from five reef systems (see Figure 1): 1) The oceanic reef systems of Ashmore Reef and 2) the Rowley Shoals; 3) the turbid and macro-tidal inshore Kimberley reef system (Adele Island, Beagle Reef and the nearshore fringing reefs within the Lalang-Garram Marine Park; 4) the fringing reefs of Gnaraloo, Quobba and Ningaloo Stations within the Ningaloo Coast World Heritage Area; and 5) Pelorus Island, midshelf central Great Barrier Reef (GBR). GBR samples were included to provide broad evolutionary and geographic context to the levels of diversity and divergence detected among reef systems in Western Australia. A total of 756 A. digitifera samples (~ 1-6 cm3) were collected from 31 sites across the four aforementioned reef systems in Western Australia (see Figure 1, Table S1), along with an additional 33 samples collected from Pelorus Island (Great Barrier Reef). Samples were identified in the field according to the morphological description provided by (Wallace, 1999). Samples were stored in 100% ethanol, subsampled and sent to Diversity Array Technology Pty Ltd. (DArT P/L) for DNA extraction, library prep, sequencing and SNP calling using the same protocol as in Thomas et al., 2020.

2. Genetic offset to climate change

Genetic offset is a term used to describe the mismatch of gene-environmental associations (GEAs) under future climate conditions (Bay, Harrigan, Underwood, et al., 2018; Fitzpatrick & Keller, 2015). This is usually characterised by the Euclidean distance between present and future biological space (Ellis et al., 2012). Under this framework, we used two model algorithms, gradient forest (GF) and generalised dissimilarity models (GDM) in the R packages ‘gradientforest’ (Ellis et al., 2012) and ‘gdm’ (Fitzpatrick et al., 2021), respectively, to describes patterns of observed genetic variation under specified climate conditions at the 26 sample sites in WA (excluding the Lalang-Garram sites). In contrast to GF which partitions the genotype data along the gradient of environmental data, GDMs are not based on machine learning and integrate distance matrices to fit gene-environmental responses using I-splines, which inform on the magnitude and slope of variables when explaining genetic turnover (Fitzpatrick & Keller, 2015; Fitzpatrick et al., 2013; Gibson et al., 2017). Once gene-environmental responses are identified at sample site locations, the models were then used to estimate regional spatial similarities in gene-environmental associations in site-neighbouring regions to predict future mismatches in these associations under climate change conditions (genetic offset) across reef systems in Western Australia, following the approach described in Fitzpatrick & Keller, 2015.

Before running ‘gradientforest’ and ‘gdm’, we identified outliers with significant GEAs using BayeScEnv (as above and excluding samples from the GBR due to high genetic dissimilarity to WA samples), which is an adapted Bayesian approach that combines FST differentiation at loci level with the selective pressure on allele frequencies driven by environmental and geomorphological conditions (de Villemereuil & Gaggiotti, 2015; Stucki et al., 2017). Loci outside the 95% false discovery rate threshold were considered outliers possibly under directional selection, and these were included in the gradient forest analysis.

Environmental variables were selected based on their importance in delineating coral growth, settlement and survival (see Table 2) (Maina et al., 2011) and can be classified into five groups; sea surface temperature (SST), SST anomalies (Zinke et al., 2018), water column optical parameters, geomorphological variables, and physical water column parameters. All variables were downscaled to the 250 m bathymetry resolution of Australia (Whiteway, 2009) using nearest neighbour resampling (Gogina & Zettler, 2010) after smoothing and completing missing environmental data using kriging interpolation (Assis et al., 2018). Once downscaled, all variables were clipped to the 0-40 m bathymetry mask, representing the zone that most photic hard corals occupy (Veron & Marsh, 1988). Prior to running BayeScEnv, variables that were correlated ³ |0.80| (Mateo et al., 2013; Senaviratna & Cooray, 2019) (Pearson correlation) with other variables at site locations were excluded to avoid overfitting, whilst retaining at least one variable from each group (see Table 2). Values of the remaining less correlated variable were extracted at each site (see Table S5) and standardised to absolute environmental distances following the BayeScEnv developers’ recommendation (Villemereuil & Gaggiotti, 2015). When extracted site variable data returned NA, values at the closest neighbouring pixel were used in further analyses. Transformed variables in association with allele frequencies of the SNP genotype data were integrated in BayeScEnv, applying default chain and model parameter settings (5,000 iterations, 20 pilot runs and 5,000 burn-in length). Posterior error probability incorporating the environmental factor (PEP g) < 0.05 was applied as recommended threshold to identify potential outlier loci or putative adaptive loci.

Putatively adaptive loci were selected for the GF analysis if they were polymorphic in more than 20% of sampled populations (Fitzpatrick & Keller, 2015) while all adaptive loci, identified using BayeScEnv are used for GDM analysis. The ‘gradientforest’ algorithm was based on 2,000 regression trees per SNP and constructed with a depth of conditional permutation adjusted to the number of variables (Fitzpatrick & Keller, 2015) and a variable correlation threshold of 0.8. GF model performance was calculated and variable importance was visualised using cumulative importance plots across individual and overall SNP’s with positive R2 value. For the GDM, the default model setting of three I-splines was used. GDM performance was assessed based on % deviance explained and the relative variable importance which was represented by the sum of I-spline coefficients (Fitzpatrick et al., 2013).

To identify the regional variation in GEA patterns and assess the future genetic offset of A. digitifera populations in WA, the study area of the 26 sites in WA was extended with a radius of 50 km {very few larvae disperse farther than 50 km (Graham et al., 2011; Jones et al., 2009; Underwood, 2009)}. The similarity in GEAs within this 50 km radius was assessed in both models based on similarities with the environmental conditions at the sampled sites (Bay, Harrigan, Underwood, et al., 2018) and visualised in Principle Component Analysis (PCA) as described in Fitzpatrick and Keller (2015). As a complementary method to determine if the spatial variable importance from the gradient forest were robust, we carried out a Sambada analysis (Sambada method and results are described in the supplementary text in the supplementary data).

Finally, we used the GF and GDM to assess the future genetic offset across the different reef systems in WA. Projected SST data from four different climate change scenarios, which we extracted from three Atmosphere-Ocean General Circulation Models (AOGCM), CCSM4, HADGEM2-ES and MIROC 5 from the CMIP 5 database (Taylor et al., 2012). We resampled future SST data to 4 km resolution using the NASA/OB.DAAC Data Analysis Software (NASA SeaDAS V 7.5.3). SST data from 2040-2050 and 2090-2100 data under RCP 2.6 (mildest scenario) and RCP 8.5 (extreme case scenario) were averaged to account for variability in future SST data. Buffer zones with a radius of 50 km of future SST data were constructed using the same downscaling and masking procedures as used for the present environmental conditions. Significant differences in genetic offset was tested between reef systems across the four climate change conditions using Kruskal-wallis non-parametric test with posthoc Dunn test with Bonferroni correction or two-way Anova with posthoc Tukey test based on the extracted Euclidean distance values within a 50 km radius around the site locations.