DOI: 10.5061/dryad.t1g1jwt4g
title: Population connectivity and genetic offset in the spawning coral Acropora digitifera in Western Australia

description: 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.
description: 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) &lt; 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.