Potential for rapid genetic adaptation to warming in a Great Barrier Reef coral

Can genetic adaptation in reef-building corals keep pace with the current rate of sea surface warming? Here we combine population genomic, biophysical modeling, and evolutionary simulations to predict future adaptation of the common coral Acropora millepora on the Great Barrier Reef. Loss of coral cover in recent decades did not yet have detectable effect on genetic diversity in our species. Genomic analysis of migration patterns closely matched the biophysical model of larval dispersal in favoring the spread of existing heat-tolerant alleles from lower to higher latitudes. Given these conditions we find that standing genetic variation could be sufficient to fuel rapid adaptation of A. millepora to warming for the next 100-200 years, although random thermal anomalies would drive increasingly severe mortality episodes. However, this adaptation will inevitably cease unless the warming is slowed down, since no realistic mutation rate could replenish adaptive genetic variation fast enough.


21
diversity in our species. Genomic analysis of migration patterns closely matched the biophysical 22 model of larval dispersal in favoring the spread of existing heat-tolerant alleles from lower to 23 higher latitudes. Given these conditions we find that standing genetic variation could be sufficient 24 to fuel rapid adaptation of A. millepora to warming for the next 100-200 years, although random 25 thermal anomalies would drive increasingly severe mortality episodes. However, this adaptation 26 will inevitably cease unless the warming is slowed down, since no realistic mutation rate could 27 replenish adaptive genetic variation fast enough.

29 30
Hot water coral bleaching, caused by global warming, is devastating coral reefs around 31 the world (1) but there is room for hope if corals can adapt to increasing temperatures (2). The 32 fact that current coral generation is suffering high mortality does not necessarily imply that the 33 next coral generation would not be better adapted. Many coral species have wide distributions 34 that span environments that differ dramatically in their thermal regimes, demonstrating that 35 efficient thermal adaptation has occurred in the past (3). But can coral adaptation keep up with 36 the unprecedentedly rapid current rate of global warming (4)? One way for corals to achieve rapid 37 thermal adaptation is through genetic rescue, involving the spread of existing heat tolerance 38 alleles from warm-adapted populations to now-warming regions via larval migration (5,6). We 39 have previously demonstrated the presence of genetic variants conferring high thermal tolerance 40 in a low-latitude A. millepora population (5). It can be hypothesized that global warming will 41 cause preferential survival of migrants from warmer to cooler locations because they will be 42 following their thermal optimum, whereas individuals migrating in the opposite direction would 43 find themselves in increasingly mismatched environments ( Fig. 1 A, B). Another likely 44 population-level effect of recent declines in coral cover (7) is a reduction in overall genetic 45 diversity, potentially limiting both the scope and the rate of adaptation. Here, we tested these 46 predictions in Acropora millepora, a common reef-building coral from the most ecologically 47 prominent and diverse coral genus in the Indo-Pacific (staghorn corals, Acropora), and used the 48 obtained demographic estimates to model the future adaptive potential of A. millepora on the 49 Great Barrier Reef (GBR).

51
Locations and genotyping 52 53 We used samples collected in 2002-2009 from five populations of A. millepora along the 54 latitudinal range of the GBR (Fig. 1 A). Environmental parameters (obtained from 55 http://eatlas.org.au/) varied widely among these locations ( Fig. 1 C). Importantly, maximum 56 summer temperature (the major cause of bleaching-related mortality) followed the latitudinal 57 gradient with one notable exception: one of the near-shore populations from the central GBR 58 (Magnetic Island) experienced summers as hot as the lowest-latitude population (Wilkie Island 59 ( Fig. 1 C).

61
We genotyped 18-28 individuals per population using 2bRAD (8) at >98% accuracy and 62 with a >95% genotyping rate. Analysis of population structure based on ~11,500 biallelic SNPs 63 agreed with previous results (9, 10) and revealed very low levels of genetic divergence, with only 64 the Keppel Islands population being potentially different from the others ( Fig. 1 D and Fig. S1).

80
Demographic subdivision and migration patterns 81 82 To more rigorously test for population subdivision and infer unidirectional migration 83 rates among populations and population sizes, we used Diffusion Approximation for 84 Demographic Inference (dadi, (11)). dadi is a method that optimizes parameters of a pre-specified 85 demographic model to maximize the likelihood of generating the observed allele frequency 86 spectrum (AFS). For two populations AFS is essentially a two-dimensional histogram of allele 87 frequencies (Fig. S2). Being a likelihood-based method, dadi can be used to compare alternative 88 models using likelihood ratio tests and Akaike Information Criterion (AIC). Most importantly, 89 unlike previously used approaches (10, 12, 13) dadi does not rely on assumptions of genetic 90 equilibrium (stability of population sizes and migration rates for thousands of generations) or 91 equality of population sizes and therefore is potentially more realistic and sensitive for natural 92 populations.

94
We used bootstrap-AIC approach to confirm that our populations are separate 95 demographic units. For each pair of populations we generated 120 bootstrapped datasets by 96 resampling genomic contigs and performed delta-AIC comparison of two demographic models, a 97 split-with-migration model and a no-split model (Fig. S3 C). The split-with-migration model 98 assumed two populations that split some time T in the past with potentially different sizes N1 and 99 N2, and exchange migrants at different rates (m12 and m21) depending on direction. The no-split 100 model allowed for ancestral population size to change but not for a population split, so the 101 experimental data were modeled as two random samples from the same population of size N. The

121
AFS-based analysis allows rigorous estimation of unidirectional migration rates between 122 populations. The classical F ST -based approach only allows estimating bi-directional migration 123 rate (13) and even that calculation has been criticized because its underlying assumptions are 124 rarely realistic (14). We determined unidirectional migration rates from the split-with-migration 125 model and estimated their confidence limits from bootstrap replicates. In theory, migration rate 126 can be confounded with population divergence time, since in the AFS higher migration often 127 looks similar to more recent divergence (15). To confirm that the model with ancient population 128 divergence and migration is preferable to the model with very recent divergence and no 129 migration, we performed the delta-AIC bootstrap comparison between these models and obtained 130 overwhelming support for the model with ancient divergence and migration (Fig. S4). Notably, 131 for all pairwise analyses migration in southward direction exceeded northward migration, and this 132 difference was significant in seven out of nine cases ( Fig. 2 A and Fig. S3A, C). Linear mixed 133 model analysis of direction dependent median migration rates with a random effect of destination 134 (to account for variation in total immigration rate) confirmed the overall significance of this 135 southward trend (P MCMC <1e-4). Full listing of parameter estimates and their bootstrap-derived 136 95% confidence limits is given in Table S1.

138
To investigate whether the southward migration bias was due to higher survival of warm-139 adapted migrants, as predicted under global warming ( Fig. 1 B

152
To determine if there were any recent changes in southward migration, we evaluated an extended 153 split-with-migration model that allowed for a change in migration over the past 75-100 years. The 154 extended model suggested some recent migration changes, including southward migration 155 increases ( Fig. S5) but once again, these changes did not correspond with an increase in migration 156 from warmer to cooler locations. We conclude that with the current data and analysis techniques 157 we cannot (yet) detect an effect of recent warming on preferential direction of coral migration on 158 the GBR.

160
Genetic diversity trends 161 162 The GBR has warmed considerably since the end of last century (18), which may have 163 already reduced genetic diversity in A. millepora populations. We used dadi to infer effective 164 population sizes, which is a measure of genetic diversity and one of the key parameters 165 determining the population's adaptive potential (19). The results of the split-with-migration 166 model ( Fig. 2 A) were consistent for all population pairs and indicated that Keppel population 167 was about one-fifth the size of others (Fig. 2 D, E). This result was not surprising since the 168 Keppel population frequently suffers high mortality due to environmental disturbances (9). We 169 also used a single-population dadi model that allowed for two consecutive growth/decline periods 170 ( (7). It must be noted, however, that recent decline in population size is hard to detect 183 using the AFS method unless the sample size is very large, since it would predominantly affect 184 the frequencies of rare alleles (15,21). Despite this shortcoming, we can conclude that genetic 185 diversity in A. millepora has not yet been strongly affected by warming over the past century, 186 although the populations appear to have been in long-term decline GBR-wide for the past several

208
We found that, with only ten thermal QTLs, under a broad range of settings for 209 heritability and plasticity the pre-adapted metapopulation was able to persist through the warming 210 for at least 20 -50 generations (100 -250 years) and, under some parameter combinations, much 211 longer ( Fig. 3 and S8). Migration substantially contributed to this persistence (

236 237
A notable tendency observed with all parameter settings was that during warming the 238 fitness (and hence the size) of adapting populations began to fluctuate following random thermal 239 anomalies, and the amplitude of these fitness fluctuations increased as the warming progressed 240 even though the amplitude of thermal anomalies did not change (Fig. 3 H)

264
Longer metapopulation persistence under low heritability and high plasticity was most 265 likely due to their enhancing effect on standing genetic variation (Fig. 5). Higher plasticity 266 promoted both higher number of variants retained in populations and larger effect sizes of these 267 variants, whereas higher heritability also led to higher number of retained variants but notably 268 smaller effect sizes (Fig. 5 A-K). It can be said that under high heritability local adaptation was 269 based on many mutations of small effect, whereas low heritability promoted adaptation involving 270 fewer mutations of larger effect. Importantly, the cumulative absolute effect of QTL variants in a 271 population was consistently higher under the setting of low heritability, despite lower number of 272 variants (Fig. 5 K, L). During warming, this variation lasted longer as the source of adaptive 273 genetic variants, enabling up to 4 o C increase in mean thermal tolerance (Fig. 3 G and S9).

274
Another important effect of higher plasticity was that it partially rescued the drop in average 275 population fitness due to low heritability (Fig. 3 B, D and S9

293 294
Potential pitfalls 295 296 There are several uncertainties in our model associated with coral biology. Below we 297 argue that, while more research is certainly needed to resolve these uncertainties, our modeling 298 was conservative overall.

300
We assumed only ten QTLs, which is likely much fewer that the actual number of 301 thermal QTLs in acroporid corals (20). Higher number of QTLs and/or their larger effect sizes 302 would promote higher genetic variation and lead to longer population persistence. We also kept 303 the distribution of QTL effect sizes narrow: with the current settings and ten QTLs, at the start of 304 simulation only about 2% of corals deviated from the mean thermal tolerance by more than 1.5 o C 305 in either direction. Such narrow variation makes adaptation to the thermal gradient of ~3 o C along 306 the GBR non-trivial, but still, at present there is no experimental data to evaluate whether even 307 such narrow variation is realistic.

309
We used effective population sizes suggested by genetic analysis as census sizes. In 310 highly fecund marine organisms census sizes tend to substantially exceed effective population 311 sizes, sometimes by orders of magnitude (29), which would strongly promote higher genetic 312 variation and population persistence. Moreover, we modeled only our five populations rather than 313 the whole GBR, which would have resulted in much higher standing genetic variation in the 314 metapopulation, promoting longer persistence.

316
As for phenotypic plasticity, in simulations shown on Figs. 3 σ = 0.5 and σ = 1 317 corresponded to 86% and 40% decline in fitness if the individual's phenotype mismatched the 318 environment by 1 o C. The existing data on coral thermal plasticity are somewhat conflicting. One 319 study shows that acroporid corals can successfully acclimatize to environments differing in 320 maximum temperatures by as much as 2 o C (30); however, another study found that coral grew 321 52-80% more slowly when transplanted among locations differing by 1.5ºC average temperature, 322 (31). Although it is not possible to directly interpret these results in terms of width of the fitness 323 function (as plasticity is encoded in our model) the former study likely supports the higher 324 plasticity setting (σ = 1) while the latter study supports σ = 0.5. It must also be noted that both 325 these studies involved in situ transplantations and hence the effect of temperature remains 326 confounded with other local fitness-affecting environmental parameters. Also, in adult corals 327 plasticity is likely lower that in larvae and recruits, which are expected to exhibit non-reversible 328 developmental plasticity associated with metamorphosis and establishment within a novel 329 environment (32). One particularly important event during this developmental transition is 330 establishment of association with local algal symbionts. Since symbionts also adapt to local 331 thermal conditions (28) this would elevate the fitness of the coral host despite possible mismatch 332 between its own genetically determined thermal optimum and local temperature, which in our

365
Conclusions 366 367 Our study provides a novel integrated empirical and modeling framework to evaluate the 368 risk of extinction in natural populations. We found that genetic diversity and migration patterns of 369 Acropora millepora were not yet strongly affected by climate change and were well positioned to 370 facilitate persistence of the GBR metapopulation for a century or more. Our results underscore 371 the pivotal role of standing genetic variation and genetic exchange in the future metapopulation 372 persistence. This implies that any intervention that would reduce this variation (for example, 373 captive breeding based on selection of only a few "winner" genotypes or clonally propagating a 374 small number of genotypes for reef restoration) is likely to have negative impact on corals' 375 adaptive potential. In contrast, efforts facilitating the spread of genetic variation, such as assisted 376 gene flow (36), could be much more helpful in the long term. With the estimated natural 377 migration rates on the order of 5-100 migrants per generation, human-assisted genotype exchange 378 could appreciably contribute to the genetic rescue without risking disruption of the natural local 379 adaptation patterns (37). However, despite good prospects for short-term adaptation, corals are 380 predicted to become increasingly more sensitive to random thermal anomalies, especially in the 381 originally warm-adapted populations. The 10-85% mortality in the Northern GBR as a result of 382 2016 bleaching event (38) could be a particularly sobering recent manifestation of this trend.

383
Finally, it is important to point out that adaptation based on genetic rescue will not save corals 384 from eventual extinction: it will only buy us some time to take action against global warming, 385 which hopefully can be stopped before corals run out of genetic variation.

389
Genotyping 390 391 This study relied predominantly on samples described by van Oppen et al (10)  ). All other samples were unrelated. We took advantage of these clonal replicates to extract 407 SNPs that were genotyped with 100% reproducibility across replicates and, in addition, appeared 408 as heterozygotes in at least two replicate pairs (script replicatesMatch.pl with hetPairs=2 option).

409
These 7,904 SNPs were used as "true" SNP dataset to train the error model to recalibrate variant 410 quality scores at the last stage of the GATK pipeline. During recalibration, we used the transition-411 transversion (Ts/Tv) ratio of 1.438 determined from the "true" SNPs to assess the number of false 412 positives at each filtering threshold (as it is expected that an increase of false positive calls would 413 decrease the Ts/Tv ratio towards unity). We chose the 95% tranche, with novel Ts/Tv = 1.451.

414
After quality filtering that restricted the calls to only bi-allelic polymorphic sites, retained only 415 loci genotyped in 95% or more of all individuals, and removed loci with the fraction of 416 heterozygotes exceeding 0.6 (possible lumped paralogs), we ended up with 25,090 SNPs. In total, 417 2bRAD tags interrogated 0.18% of the genome. The genotyping accuracy was assessed based on 418 the match between genotyped replicates using script repMatchStats.pl. Overall agreement 419 between replicates was 98.7% or better with the heterozygote discovery rate (fraction of matching 420 heterozygote calls among replicates) exceeding 96%.

422
Genome-wide genetic divergence 423 424 To begin to characterize genome-wide divergence between populations we used pairwise 425 genome-wide Weir and Cockerham's F ST calculated by vcftools (45)

487
Coral-specific biological parameters for A. millipora included relative adult density 488 (dependent on the habitat), reproductive output, larval spawning time and periodicity (e.g.,

489
Magnetic Island populations spawn a month earlier than the other GBR sites (52)), maximum 490 dispersal duration, pre-competency and competency periods, and larval mortality (53, 54). The 491 spatially explicit dispersal simulations model the dispersal kernel (2-D surface) as a 'cloud' of 492 larvae, allowing it to be concentrated and/or dispersed as defined by the biophysical parameters.

493
An advection transport algorithm is used for moving larvae within the flow fields (55).

495
Simulations were carried out by releasing a cloud of larvae into the model seascape at all 496 individual coral reef habitat patches and allowing the larvae to be transported by the currents.

497
Ocean current velocities, turbulent diffusion, and larval behavior move the larvae through the 498 seascape at each time-step. Larval competency, behavior, density, and mortality determine when 499 and what proportion of larvae settle in habitat cells at each time step. When larvae encounter 500 habitat, the concentration of larvae settling with the habitat is recorded at that time-step. From the 501 dispersal data, we derived the coral migration matrix representing the proportion of settlers to 502 each destination patch that came from a source patch, which is analogous to the source 503 distribution matrix (56) and is equivalent to migration matrices derived from population genetic 504 analysis. It is important to note that migration matrices extracted for the field sites represent the 505 potential migration through all possible stepping-stones.

507
Metapopulation adaptation model 508 509 The model was implemented in SLiM (22) we modeled only ten QTLs with normal distribution of effect sizes with a standard deviation 534 of 0.2 o C. With ten QTLs, this setting implied that at the start of simulation only about 2% of 535 corals deviated from mean thermal tolerance by more than 1.5 o C in either direction. Since 536 thermal differences between our populations exceeded 3 o C, this narrow variation made local 537 adaptation rather non-trivial.

538
-Dominance of QTLs (set to 0.5 in our simulation).

539
-Phenotypic plasticity: standard deviation of the Gaussian curve describing decline in fitness 540 away from phenotypic optimum. We modeled three plasticity settings, 0.5, 1 and 2, which 541 corresponded to 86%, 40% and 13% fitness drop when the environment mismatched 542 phenotypic optimum by 1 o C.

543
-Non-heritable phenotypic component: standard deviation of a normal distribution with mean 544 zero from which a random value is drawn to be added to the sum of QTL effects when 545 computing phenotype. Setting this parameter to zero corresponds to trait heritability of one.

546
Higher values of this parameter imply heritability less than one; however, the exact value of 547 heritability (the proportion of phenotypic variation explained by genetics) could still vary 548 depending on the extent of genetic variation. 549 -Mutation rate. It was either set to zero to explore the role of standing genetic variation or 550 varied between 1e-6 and 1e-4 per QTL per gamete. This range covers and exceeds the range 551 of trait-level deleterious mutation rates observed in humans (57). Therefore, values at higher 552 end of the range most likely strongly over-estimate the rate of adaptive mutations, which was 553 done deliberately to show that no realistic mutation rate could help sustain genetic variation 554 in the face of strong selection by warming.

571
All combinations of parameter settings were run ten times to ensure consistency. We 572 found that with population sizes in thousands, such as in our case, the results were very 573 consistent among independent runs. We therefore did not aggregate results over many replicated 574 runs but show one randomly chosen run for each tested parameter combination. Esd implies lower heritability) and phenotypic plasticity (σ, standard deviation of the Gaussian slope of 818 fitness decline away from the phenotypic optimum, in degrees C). 819 820