Characterizing the adaptive potential of north-western Pacific reef corals to empower conservation strategies: a role for seascape genomics

Coral reefs are suffering a major decline due to the environmental constraints imposed by climate change. Over the last 20 years, three major coral bleaching events occurred in concomitance of anomalous heat waves, provoking a severe loss of coral cover world-wide. Recent works, however, reported of reefs at which recurrent bleaching events resulted in increased resistance to thermal stress, suggesting that adaptation might have occurred. The conservation strategies for preserving the reefs, as they are conceived now, can not cope with global climatic shifts. In this regard, researchers advocated the set-up of a preservation framework to reinforce coral adaptive potential. The main obstacle to this implementation is that studies on coral adaption are usually hard to generalize at the scale of a reef system. In this work, we ran a seascape genomics and connectivity analysis to characterize the adaptative potential of a flagship coral species of the Ryukyu Archipelago (Japan). By associating genotype frequencies with descriptors of historical environmental conditions, we discovered six genetic polymorphisms that might promote resistance against thermal stress and predicted their presence across the whole study area. Furthermore, we confronted genetic distances between reefs with their spatial separation according to sea currents to predict connectivity patterns across the region. The cross-talk between the results of these analyses portrayed the adaptive potential of this population: we were able to identify reefs carrying potential adaptive traits and how do they disperse to neighbouring reefs. This information was summarized by objective, quantifiable and mappable indexes covering the whole region that can be extremely useful for prioritization in conservation plans. This framework is transferable to any coral species and on any reef system and therefore represents a valuable tool for empowering preservation efforts dedicated to the protection of coral reef in warming oceans.


Introduction
Coral reefs are suffering a severe decline due to the effects of climate change (Hughes et al., 2017). Loss of reef is already showing catastrophic consequences for marine wildlife depending on these structures (Pratchett, Thompson, Hoey, Cowman, & Wilson, 2018) and disastrous aftermaths are expected for human economies as well (Moberg & Folke, 1999). One of the major threat to the persistence of this ecosystem is coral bleaching (Bellwood, Hughes, Folke, & Nyström, 2004): a physiological response induced by environmental stress that provokes hard-skeleton coral, the cornerstone of the reef, to separate from the symbiotic microbial algae living within their cells (Mydlarz, McGinty, & Harvell, 2010). These algae, called Symbiodiniacea (LaJeunesse et al., 2018), are crucial for the host sustainment and therefore a persistent bleached state leads coral to death (Mydlarz et al., 2010). Three mass bleaching events occurred over the last 20 years in concomitance with anomalously high sea water temperatures, reason why thermal stress is now considered the main driver of this phenomenon (Hughes et al., 2017). These events struck world-wide and lead to local losses of coral cover up to 50% (Hughes et al., 2018(Hughes et al., , 2017. Besides, coral bleaching was found to be associated with other stressors linked to human activity, such as ocean acidification, water eutrophication, sedimentation and overfishing ( Conservation efforts facing the plague of coral bleaching aim to restore reefs that underwent severe losses and to limit the impact of future damages (Baums, 2008;Bellwood et al., 2004;Young, Schopmeyer, & Lirman, 2012). For these purposes, two tools are nowadays usually used: the establishment of marine protected areas (MPAs) and the set-up of coral nurseries (Baums, 2008;Bellwood et al., 2004;Young et al., 2012). MPAs are designed zones in which human access and activities are severely restricted in order to alleviate the effects of local anthropogenic stressors (Lester et al., 2009). Coral nurseries are underwater gardens of transplanted colonies that can then be employed to restore damaged reefs (Baums, 2008;Young et al., 2012). In both cases, researchers advocated the use of design systems accounting for demographic connectivity so that the position of a conservation measure can promote resistance and resilience even on neighbouring sites (Baums, 2008 To this end, seascape genomics tools are likely to play an important role. Seascape genomics is the marine counterpart of Landscape genomics, a branch of population genomics that investigates adaptive potential in field-based experiments (Balkenhol et al., 2017;Rellstab, Gugerli, Eckert, Hancock, & Holderegger, 2015). Samples are collected across a landscape and then genotyped using next-generation-sequencing techniques, providing thousands of genetic variants (Rellstab et al., 2015). Simultaneously, the environmental characterization of the study area is carried out mainly based on remote-sensing data describing specific local climatic conditions (Rellstab et al., 2015). Genomics and environmental information are then combined in order to detect genetic polymorphisms associated with peculiar conditions (i.e., potentially adaptive traits; Rellstab et al., 2015). This approach has already been applied to many land species and is increasingly used to analyse marine animals in what is referred to as seascape genomics (exhaustively reviewed in Riginos, Crandall, Liggins, Bongaerts, & Treml, 2016).
As far as we know, no seascape genomics experiment was applied to reef corals yet. In fact, adaptation of these species has been mostly studied by performing transplantation assays followed by conditioning in aquarium, a time and resource-demanding approach that often focuses on a couple of contrasting reef ( (Balkenhol et al., 2017). Consequently, the main (adaptive) and collateral (connectivity) information generated by this method provides a valuable contribution for preservation strategies (Baums, 2008;Krueck et al., 2017;Palumbi, 2003).
In this work, we applied a seascape genomics framework to detect coral reefs carrying potential adaptive traits to show how conservation policies could implement such findings. Our study focuses on Acropora digitifera of the Ryukyu Archipelago in Japan (Fig. 1), an emblematic species of the Indo-Pacific and flagship organism for studies on corals genomics (Shinzato et al., 2011). We first analysed the convergence between genomic and environmental information i) to detect loci potentially conferring a selective advantage and ii) to develop a model describing connectivity patterns. Next, we took advantage of these findings to indicate which reefs were more likely to be carrying adaptive traits and to evaluate how well such sites are interconnected with the rest of the reef system. Finally, we propose an approach to implement the results obtained in conservation planning. In general, our work provides tools at the interface of conservation genomics and environmental sciences likely to empower preservation strategies for coral reefs.

Materials and methods
The input genomics data used in this paper comes from a published work from Shinzato et al. (2015) describing the population structure of A. digitifera of the Ryuyku Archipelago. Our framework is structured in two axes of analysis and prediction, one focusing on the presence of adaptive traits (seascape genomics) and the other on connectivity (Fig. 2). The seascape genomics analysis ( Fig. 2A) crosses genomics data with environmental information at sampling sites to discover potentially adaptive genotypes. The model describing this relationship is then used to predict, at the scale of the whole study area, the probability of presence of adaptive genotypes (Fig. 2B). In the connectivity study (Fig. 2C), we built a model describing how sea-distances between sites (computed out of remote sensing data) correspond to genetic-distances. This model is then used to predict connectivity at scale of the whole study area (Fig. 2D). Finally, the predictions on where the adaptive traits are more likely to appear and on how the reef system is interconnected allow to assess the adaptive potential across the whole study area (Fig. 2E).

Genomic dataset
The genomic data used come from a publicly available dataset consisting of 155 geo-referenced colonies of A. digitifera from 12 sampling locations (13±5 colonies per site) of the Ryukyu Archipelago in Japan ( Fig. 1 and raw sequencing reads were aligned Figure 2. Workflow of the methods. The starting point for the analysis is the generation of genetic data describing the genotypes observed at different sampling locations (in this example, 4 sampling sites). In the seascape genomics analysis (A) this data is crossed with environmental information to discover genotypes whose frequencies associate with specific climatic conditions (ENV). Such genotypes are defined a potentially adaptive. The model describing this link is then fit with environmental data at the scale of the whole reef system (B), to predict the probability of presence of the adaptive genotype (green: high probability, red: low probability). The genetic data is also crossed with seacurrent information to build a connectivity model (C) describing how seadistances between sampling sites correspond to genetic-distances. This model is fit with sea-distances between all the reefs of the study area to predict (D) patterns of connectivity from (outbound) and to (inbound) each reef (green: high connectivity, red: low connectivity). Finally, predictions on the presence of adaptive genotypes and connectivity patterns are combined to assess the adaptive potential across the study area (E): reefs that are connected with sites likely to carry the adaptive traits will be have a high adaptive potential (green), vice versa for those isolated (red).

Genetic Data
Probability This pipeline produced the filtered genotype matrix consisting of 136 individuals and 15,516 SNPs.

Environmental data
Eleven georeferenced datasets describing atmospheric and seawater conditions were retrieved from publicly available resources (EU Copernicus Marine Service, 2017; NASA, 2016; National Oceanic and Atmospheric Administration, 2017; Tab S1). These remote sensing derived datasets provided monthly or daily variables measured across the whole study area, with a spatial resolution ranging from 25 km to 4 km (Tab S1). We processed these variables in the R environment using the raster package (v. were assigned to each reef-cell using the QGIS zonal statistics tool.

Seascape genomics
We performed the genotype-environment association analysis using the logistic regression  Jombart, 2008). This procedure highlighted a main separation between two groups of samples along the first discriminant function (Fig. S1), which was therefore used as co-variable in association models. We then ran the pipeline for the discovery of adaptive signals described hereafter and summarized in a schematic view in

Probability of Presence of Adaptive Genotypes
The results of the seascape genomics analysis were then used as starting point for predicting the probability of presence of adaptive genotypes at the scale of the whole Ryukyu Archipelago (Fig. 2B). These genotypes can be theoretically associated to any environmental descriptor but, since most of the putative adaptive genotypes were found to be associated with DHW frequency (Tab. 1), we designated this variable as the environmental pressure of interest.

Sea current data
Daily records of sea surface current were retrieved from publicly available databases (zonal and meridional surface velocities from the global-reanalysis-phy-001-030 product; EU Copernicus Marine Service, 2017) and used to compute the direction and speed of currents in the R environment using the raster library. These day-by-day calculations were then piled-up to retrieve, for each pixel of the study area (~0.083°), the cumulated speed (the conductance) toward each of the eight neighbouring pixels. This information was used to calculate dispersal costs (the inverse of conductance) and was summarized in a transition matrix in the format of the R gdistance package (v. 1.2, van Etten, 2018). For the connectivity analysis (Fig. 2C), the transition matrix was used to calculated sea-distances (i.e. the least-cost path) between sampling sites of the genotyped colonies. For the connectivity predictions (Fig. 2D), we calculated the sea-distances between all the reefs of the study area, whose coordinates were retrieved using the same approach previously described in the environmental variables section.
Importantly, for every couple of reefs (for instance reef1 and reef2) two sea-distances were computed, one per direction (i.e. from reef1 to reef2 and from reef2 to reef1).

Connectivity analysis
Genomics data was processed by the means of a principal component analysis (PCA; Fig. S3) on the filtered genotype matrix: the values of the PCs were used to calculate the Euclidean distances between each couple of samples that were subsequently averaged by site in order to obtain GDs between reefs. Next, we built a linear model (hereafter referred to as the connectivity model) to estimate GDs out of log-transformed sea-distances and ran a t-test to check for the statistical significance of the relationship. Since we wanted to consider the genetic distances within reef (i.e. sea-distance = 0) we added +0.005 to all the sea-distances prior to the log-transformation to avoid obtaining non-definite values (i.e. the logarithm of zero). As a comparison, we also built a connectivity model using log-transformed aerial distances (i.e. Euclidean distance of coordinates) as independent variable while maintaining GD as response variable. We then confronted models based on sea-and aerial-distance by comparing the coefficients of determination (R2) and the Akaike information criterion (AIC; Bozdogan, 1987).

Connectivity Predictions
We then used the connectivity model to predict GDs between any couple of reefs of Ryukyu Archipelago. This was done by translating sea-distances in GDs using the connectivity model trained during the connectivity analysis. Importantly, since there were two sea-distances for every couple of reefs (i.e. from reef1 to reef2 and vice versa), two GDs were computed as well.
Based on these predictions, we were able to calculate, for each reef-cell, two indices (Fig. S4): -outbound connectivity (OC; Fig. S4A): the % of cells reachable from the focal cell within a GD-threshold (tGD).
-inbound connectivity (IC; Fig. S4B): % of cells that can reach the focal cell within a tGD.
Since these indices were conceived to display contrasts in connectivity among reefs of the Ryukyu Archipelago, we looked for the tGD value maximizing their variances (Fig. S5). The values corresponding to the 25%-quantile and the 60%-quantile of the GD distribution were those boosting the differences in both connectivity indices. To avoid the artificial bias of border effects, we preferred working at a more local scale and set tGD at the 25%-quantile of the GD distribution (i.e. tGD=79.861).

Evaluation of the adaptive potential
The adaptive potential was evaluated by combining the predictions on the presence of adaptive traits and on connectivity patterns across the reefs of Ryuyku Archipelago (Fig. 2E). The previously computed average probability of carrying adaptive traits (Pa) is used to identify the putatively adaptive reefs. These are the reefs approaching the maximal values of Pa (Pa>0.6).
We then used the GDs calculated during the connectivity analysis to retrieve, for each reefcell, the GD from the closest putatively adapted reef. This results in an index describing the adaptive potential of every reef of the study area: the distance from the closest adaptive reef.

Seascape Genomics
We detected 15 significant genotype-environment association-models (A1-A15; qG and qW<0.001; Tab.1), referring to genetic markers from 6 distinct scaffolds of the A. digitifera reference genome. Among them, one association was related to SST variation in April, one to changes in sea current velocity in May and 13 to DHW frequency (or to a highly correlated variable such as variation in cloud fraction in August).
The SNP in the significant association-model A1 is related to SST standard deviation in April and lies within a genomic window (± 250 kb) counting 18 predicted genes, 7 of which are annotated with a known function (Box S1). Association-models A2 to A5 include four SNP  S6). Association-models A8 to A14 refer to seven SNPs present in the same scaffold (~195 kb long). This region contains five putative genes of which 2 annotated by similarity (Box S6). In association-model A15, finally, the SNP associated with DHW frequency is located on a ~75 kb scaffold that contains only one putative gene predicted as ATP-dependent DNA helicase RecQ-like (LOC107328703; Box S7). The seascape genomics analysis using the SamBada method detected 15 significant association-models (qG and qW<0.001). This table shows, for each association-model, the genomic position (scaffold of the reference genome and relative situation), the genotype involved (0 and 2=homozygous, 1=heterozygous), the concerned environmental variables (in brackets is displayed a variable from the same group of correlated descriptors that resulted in a model with a higher G-score, if any), the q-values associated with the G-score (qG) and the Waldscore (qW), the number of predicted genes in a ± 250 kb window around the SNP (with the number of SNPs annotated by similarity with a known gene, in brackets) and the annotations from the reference genome that are in the neighbourhood (± 250 kb) of the SNP and that corroborate an adaptive role against the respective environmental variable.

Probability of Presence of Adaptive Genotypes
The significant association-models of the seascape genomics analysis (Tab. 1) were used to transform DHW frequency values in probabilities of presence of adaptive genotype (Fig. S7).
The average of these probabilities (Pa) ranges from 0.40 to 0.62 (Fig. S7). Potentially adapted reefs (i.e. reefs with Pa>0.6) are located in the islands of Okinawa and Miyako, in the remote Japanese islands east of the Ryukyu Archipelago and in the Northern part of Philippines (Fig.   S7).

Connectivity Analysis
The two dispersal models respectively based on sea-and aerial-distances resulted in a significant association between genetic-distance (GD) and spatial separation (p<10 -9 , Fig. S6).
We found that the model based on sea-distances better explained the variation in GD (R2=0.57) than the aerial-distances model (R2=0.45). The AIC is lower for the model based on seadistances (AIC=315) than for the aerial-distances model (AIC=343).

Connectivity Predictions
The reefs around the islands of Tokara were those with the highest values of IC (0.49±0.01; Evaluation of the adaptive potential

Adaptation to thermal stress
Thermal stress is expected to be one of the major threat to coral reef survival and the quest for  The map displays the predicted genetic-distance (GD) from the closest reef expected to carry adaptive genotypes against thermal stress. Adaptive sites are designed as those having a cumulative probability of carrying adaptive genotypes (Pa, Fig. S7) higher than 0.6. The GDs are computed out of sea-distances using the dispersal models developed in the connectivity analysis.

Miyako
Okinawa Amami Tokara Osumi on the expression of neighbouring genes are more frequent and can influence loci at 1-2 Mb of distance (Brodie et al., 2016). The fragmentation of the reference genome forced us to limit our search window to ±250 Kb around each SNP and still we found several annotations corroborating a response to heat stress.
The SNP in association-model A6 is related to "JNK1/MAPK8-associated membrane protein", the closest gene with a predicted function (Box S5). Mitogen-activated-protein kinases Seascape/Landscape genomics studies are susceptible to high false discovery rates, especially when neutral genetic structure is not accounted for (Rellstab et al., 2015). To take this element into account, we processed data following a conservative pipeline and used models explicitly integrating demographic processes. We also set a restrictive threshold to filter significant association-models (i.e. q<0.001 in two distinct tests). Most of the significant associationmodels found were linked to DHW, most probably because this variable represents one of the main constraint to coral survival (Hughes et al., 2017). It is also possible that the initial application of a sampling scheme adapted to seascape genomics (unlike the one used by

Patterns of connectivity
Understanding connectivity patterns across a reef system is essential in order to plan effective 1 management strategies. This is usually carried out through the quantification of the dispersal 2 trajectories of larvae (Storlazzi, van Ormondt, Chen, & Elias, 2017) or through the study of the 3 population structure using genetic data (e.g. Lukoschek et al., 2016). The first approach has the 4 main drawback of being difficult and expensive because it requires an a priori knowledge of 5 the dispersal of the studied species (McCook et al., 2009). In population genetics data, on the 6 other hand, sampling often occurs at a few locations and the results are therefore hard to 7 generalize to a whole seascape extent (Beger et al., 2014). The method we applied combines 8 these two approaches. Firstly, genetic data are used to calculate genetic distances between the 9 sampled locations. The latter were then used to train dispersal models, which could in turn 10 predict genetic distances at the scale of the whole reef system. 11 In model training we found that sea-distances clearly represented a better indicator of spatial 12 separation, compared with aerial distances (Fig. S6). This happens probably because coral 13 dispersal is driven by the water flow (Paris-Limouzy, 2011), which is strongly asymmetrical 14 in this region (north-east oriented) because of the Kuroshio current (Nishikawa, 2008). By 15 fitting this dispersal model with sea-distances among any reef of the region we were able to 16 evaluate the patterns of connectivity. Reefs in the extreme North (Osumi) face both inbound 17 and outbound isolation (Fig. 3), a situation that might lead to inbreeding depression (Keller & 18 Waller, 2002). Reefs in the extreme south-west of the system (Yaeyama) appears also 19 endangered: despite their proximity to Miyako, the strong north-eastward directionality of sea-20 currents contrasts the arrival of recruits from the rest of the Archipelago. 21 The dispersal framework we propose here is potentially applicable to any coral species for 22 which geo-referenced genetic data are available and the indices associated can be extremely 23 valuable in management planning (Beger et al., 2014). However, a substantial part of GD a few loci), which is less informative for adaptive processes, compared to genomic data. Here 39 we showed how a combination seascape genomics and connectivity analysis makes it possible 40 to assess the adaptive potential at the scale of a reef system. 41 In the specific case of the Ryukyu Archipelago, we found that reef from Miyako and Okinawa 42 islands are those more likely to carry adaptive traits against thermal stress (Fig. 4) thermal stress before and therefore are not predicted to show the same adaptive capacities 47 (Fig. 4). 48 The indicators on adaptive potential and connectivity produced in this work provide valuable 49 descriptors of the state of the reef. Furthermore, they can be used for predicting the expected 50 impact of distinct preservation measures and support therefore spatial planning in 51 conservation (Box. 1). It is unknown whether the connectivity patterns and the adaptive signals 52 observed in A. digitifera are representative for those of other taxa. In this concern, we advocate 53 that future research should focus on a multi-species implementation of the methods proposed 54 here. 55

Box. 1 Adaptive potential in conservation scenarios.
In this box we show how information on adaptive potential can be used to compare different conservation scenarios. We focus on two typical preservation initiatives, marine protected areas (MPAs) and coral nurseries (CNs). For each of them, we compare how the set-up of the measure at two different locations differentially drives the adaptive potential of the reef system. The starting point for the predictions is the map showing the distribution of the probability of presence of adaptive genotypes (Fig. S7). We modify the values in this map to imitate the effects of preservation strategies. In the case of MPAs, we maintain potentially adaptive reefs (Pa>0.6) only at locations where the protected area is established. The underlying idea is to emulate the way MPAs reduces coral loss due to local stressors (Selig & Bruno, 2010). In the CNs scenarios, potentially adaptive reefs (Pa>0.6) are added where absent. Here we follow the good practices for this kind of measure that suggest to transplant adapted colonies at sites lacking them (Baums, 2008). After modifying the values of Pa, we re-calculate the genetic distance (GD) from adaptive sites (Fig. 4). This is the variable that we use to describe adaptive potential. We then compare differences in this value (ΔGD) between different conservation scenarios to investigate beneficial and negative effects.

Marine Protected Areas
We compare two MPAs established in Miyako-Yaeyama area (MPA1) or in Okinawa (MPA2), both covering the same surface (~15'000 km 2 ). The overall difference in GD (ΔGD=+0.17) indicates that MPA2 globally reduces GDs from potentially adapted reefs more than MPA1. This is probably due to the fact that MPA2 includes more reefs that are putatively adaptive (Fig. S7) and is better connected to the reef-dense area around Amami. MPA1, on the other hand, is far from this reef-dense area and shows beneficial effects only in Yaeyama.

Coral Nurseries
Two CNs emplacements were compared: in Amami (CN1) or in Tokara (CN2). The overall difference in GD (ΔGD=-0.28) suggests a more beneficial effect on adaptive potential under CN2, rather than CN1. In fact, the existence of a coral nursery in Amami (CN1) could provide adaptive genetic material to a larger number of reefs, compared with a set-up established in Tokara (CN2), because of a higher outboundconnectivity potential (Fig. 3b).

Conclusions 57
This study highlights the value of the seascape genomics approach to support the conservation 58 of corals. We applied this approach to a flagship coral species of Ryukyu Archipelago and 59 identified genetic variants that might underpin adaptation to thermal stress. Coupling this 60 information with a genetic analysis of connectivity enabled the evaluation of the adaptive 61 potential at the scale of the whole study area. The outputs of this analysis are quantitative 62 indices likely to support objective prioritization of reefs in conservation plans. This framework 63 is transferable to any coral species on any seascape and therefore constitutes a useful 64 conservation tool to take into account and evaluate the genomic adaptive potential of coral 65 reefs worldwide. 66 67 68