Genomic models predict successful coral adaptation if future ocean warming rates are reduced

See allHide authors and affiliations

Science Advances  01 Nov 2017:
Vol. 3, no. 11, e1701413
DOI: 10.1126/sciadv.1701413


Population genomic surveys suggest that climate-associated genetic variation occurs widely across species, but whether it is sufficient to allow population persistence via evolutionary adaptation has seldom been quantified. To ask whether rapid adaptation in reef-building corals can keep pace with future ocean warming, we measured genetic variation at predicted warm-adapted loci and simulated future evolution and persistence in a high-latitude population of corals from Rarotonga, Cook Islands. Alleles associated with thermal tolerance were present but at low frequencies in this cooler, southerly locality. Simulations based on predicted ocean warming in Rarotonga showed rapid evolution of heat tolerance resulting in population persistence under mild warming scenarios consistent with low CO2 emission plans, RCP2.6 and RCP4.5. Under more severe scenarios, RCP6.0 and RCP8.5, adaptation was not rapid enough to prevent extinction. Population adaptation was faster for models based on smaller numbers of additive loci that determine thermal tolerance and for higher population growth rates. Finally, accelerated migration via transplantation of thermally tolerant individuals (1 to 5%/year) sped adaptation. These results show that cool-water corals can adapt to warmer oceans but only under mild scenarios resulting from international emissions controls. Incorporation of genomic data into models of species response to climate change offers a promising method for estimating future adaptive processes.


Predicting future species response to climate change requires detailed knowledge of the link between organismal success and environment. Under most predicted scenarios, local temperatures will surpass current thermal tolerance limits of many species (1). However, many traditional climate envelope models for predicting species responses to climate change assume that thermal limits are static over time (2). Although more complex models have attempted to integrate variation in tolerance in the form of plasticity or evolutionary adaptation into estimates of vulnerability to climate change (3, 4), these models are often hindered by a dearth of information about the genetic underpinnings of climate-related traits, which often define the rates and limits of adaptation.

Increasing observations suggest that genetic variation associated with thermal tolerance could provide the raw material necessary for adaptation to increasing temperatures (5, 6). Population genomic studies from a wide array of taxa have consistently found allelic variation associated with climate traits (79). These genetic variants tend to be numerous and distributed across a broad array of gene functions, suggesting that thermal tolerance may be a cumulative product of many small-effect loci (10, 11). Because thermal tolerance is often not determined by one or a few loci of large effect and thus is unlikely to be fixed between populations, most populations likely contain some standing variation in genetic thermal tolerance. As climate change proceeds, natural selection might increase the frequency of these thermal tolerance genotypes through differential survival and propagation of preadapted genomes. The main questions that revolve around this process are whether the rates and limits of adaptation are high enough to (i) assure a population’s survival, (ii) delay its environmental demise, or (iii) neither.

In the ocean, reef-building corals are among the organisms most vulnerable to rising temperatures. Although corals inhabit a wide variety of environments, they are thought to live at just 1° to 2°C below their upper thermal limits (12). Despite this, a number of studies have provided evidence for adaptation of coral populations to local temperatures (13, 14) and for heritable variation in thermal tolerance in corals (11, 15). We capitalize on a previous study (11) on the coral Acropora hyacinthus in which we identified alleles at many loci associated with thermal tolerance across a temperature gradient in American Samoa. Here, we focused on a high-latitude population of the same species in Rarotonga, Cook Islands, asking whether this population harbored the same warm-adapted alleles found in our previous study, thus providing the potential for adaptation of this population to warmer temperatures. Next, we used a coupled environmental, evolutionary, and demographic simulation to predict survival or extirpation under future climate change scenarios. Finally, we tested whether migration of heat-tolerant genotypes, for example, through outplanting, could increase the likelihood of population persistence. Our results are preliminary because our knowledge of genotype-phenotype effects and population dynamics in this system is imperfect. Nevertheless, we develop a general framework by which questions of future evolutionary rate can be answered and show that evolution in this population can substantially increase future survival.


Presence of heat-adapted alleles in a high-latitude population

We previously identified 114 alleles associated with thermal tolerance in A. hyacinthus on Ofu Island, American Samoa, and compared allele frequencies of these corals from two pools, less than 1 km apart, with differing temperature profiles (11, 16). The corals in these pools represent the same broad gene pool; but they are affected by strong heat-related natural selection and diverge at these 114 loci (11, 17). Here, we focus on a high-latitude population of the same species in one of the most southerly coral reefs in the Pacific, Rarotonga, Cook Islands. Large A. hyacinthus colonies are sparsely distributed across the fringing-reef and back-reef areas in Rarotonga. The island is more than 200 km from even the smallest nearby island and thus likely experiences restricted dispersal across ocean gaps. As a result, coral migration in and out of Rarotonga may be limited, and we treat this population as effectively ecologically closed.

First, we tested whether this isolated, high-latitude population harbored the same warm-adapted alleles found in the HV (high temperature variation) pool of Ofu, potentially allowing for adaptation. Combining RNA sequencing (RNA-seq) data sets from Ofu and Rarotonga, we were able to examine 20,883 single-nucleotide polymorphisms (SNPs). Of these, 6451 SNPs occurred at a low frequency (<5%) (Fig. 1B). Overall, among SNPs at a moderate frequency (>5%), fixation index (FST) was higher between Ofu and Rarotonga than between the HV and MV (moderate temperature variation) pools on Ofu Island (Fig. 1A). These differences across archipelagos were also reflected in a principal components analysis showing that the largest proportion of the variation (PC1, 4.9%) corresponded to the differences between Ofu and Rarotonga (Fig. 1C).

Fig. 1 Patterns of polymorphisms across populations.

(A) FST values between Ofu Island, American Samoa and Rarotonga, Cook Islands show a higher frequency of large FST than comparisons within Ofu Island. Inset magnifies the tail of distribution. (B) The Rarotonga population has a greater number of rare alleles (<5% frequency overall) than Ofu. (C and D) Patterns of differentiation are also reflected in a principal components analysis of SNPs from A. hyacinthus colonies from Rarotonga, Cook Islands (Rar) and two sites in Ofu Island, American Samoa: an HV pool and an MV pool. The first principal component (PC1) across all SNPs divides groups by population (C), whereas the first principal component across candidate SNPs for heat tolerance divides groups by microclimate (D).

However, when we examine only 114 SNPs that were associated with thermal environment in our previous study (11), we see a different pattern. We had sufficient coverage in Rarotonga to examine 65 of 114 thermal tolerance SNPs (table S1): Across these, frequencies were more similar between the cooler MV pool and Rarotonga than between the warmer HV pool and Rarotonga, although all locations were significantly different from one another (Fig. 1D; P < 0.01, Tukey post hoc test). In the HV pool, the fraction of putative thermal tolerance alleles (Tprop) for each individual coral across 114 SNPs ranged from 0.4 to 0.58 (mean, 0.49). Tprop was lower in the MV pool (0.11 to 0.28; mean, 0.18) and intermediate in Rarotonga (0.19 to 0.38; mean, 0.28). We would not expect allele frequencies in Rarotonga to be higher than those in the MV pool on the basis of thermal environment alone, so it is likely that other facets of selection and drift act on these loci. Nevertheless, the MV pool and Rarotonga both show much lower allele frequencies at these SNPs than the HV pool.

Probability of population persistence under climate change

We conducted population genetic simulations to estimate evolution at putatively adaptive loci using a conversion from genotype to temperature tolerance based on the differences between the HV and MV pools in Ofu (Fig. 2; see Materials and Methods). We used ocean warming conditions predicted for four representative concentration pathways (RCPs) developed for the Intergovernmental Panel on Climate Change (IPCC) using remote sensing data from Rarotonga and population growth rates estimated from a previous study (18), which measured recruitment on an isolated reef after disturbance.

Fig. 2 Genomic prediction of Rarotonga population trajectories under RCP emissions scenarios.

(A) Temperature of the hottest month from NOAA GFDL ESM2M (National Oceanic and Atmospheric Administration Geophysical Fluid Dynamics Laboratory Earth System Model 2M) SST predictions for Rarotonga. (B) Population size trajectories. For each RCP, 10 simulated populations are shown in transparent colors, and the mean population size is shown in bold. (C) Frequency distributions for thermal tolerance SNPs identified by Bay and Palumbi (11). Black line reflects initial frequency distribution, and dark and light blue lines show RCP2.6 and RCP4.5, respectively. The other RCP scenarios are not shown because the simulated populations went extinct.

Populations simulated under RCP2.6 and RCP4.5 persisted through 2100 (Fig. 2). In RCP2.6, population sizes tracked changes in temperature, even declining when temperatures showed strong annual swings (for example, after 2075). The fraction of thermal tolerance alleles (Tprop) also tracked the environment, increasing from an initial mean of 0.28 to 0.35 under RCP2.6. In RCP4.5, populations did not initially track quickly warming conditions, falling below 600 individuals between 2035 and 2065. Strong selection led to increased adaptation (Tprop changed from 0.28 to 0.43) and resulted in population growth. These adjustments in allele frequencies and population sizes suggest that the rate of adaptation largely kept pace with warming rates of RCP2.6 and RCP4.5. By contrast, under RCP6.0, 6 of 10 simulated populations went extinct, and the remaining populations had less than 50 individuals by the year 2100. Under RCP8.5, all simulated populations went extinct. Under RCP6.0, there is an initial decline and recovery similar to that of RCP4.5, but as temperatures climbed above 29°C, the populations collapsed. Under RCP8.5, populations also declined precipitously after temperatures reached an average of 29°C by about 2040. These results from our probabilistic model were nearly identical to those using a second set of simulations from an individual-based model (fig. S1). In addition, our results were robust to the choice of climate model; using the mean sea surface temperature (SST) output across all 35 IPCC class models from the Fifth Assessment Report (AR5) yielded qualitatively similar results, although shapes of the curves differed slightly (fig. S2).

Recent efforts in coral conservation have included a variety of outplanting practices, with suggestions that moving corals from high to low tolerance locations could enhance local persistence. Because the efficacy of assisted gene flow in climate adaptation has seldom been tested, we conducted simulations using RCP6.0 to determine whether migration of heat-tolerant genotypes could increase the likelihood of population persistence. We found that yearly transplants of 10 heat-tolerant individuals, with allelic variation drawn from the HV pool, increased the population size at 2100 from near extinction (<25 on average) without migration to a stable and growing average of 399 individuals (Fig. 3). Larger migration rates had more marked effects: Adding 25 and 50 individuals/year resulted in average 2100 population sizes of 600 and 728, respectively. This was not a purely demographic effect—addition of non–heat-tolerant individuals (with Rarotonga allele frequencies) still led to extinction (fig. S3).

Fig. 3 Assisted gene flow affects adaptation to climate change.

(A and B) Four different levels of migration (0, 10, 25, and 50 migrants/year) simulated under RCP6.0. Ten simulated trajectories (A) and associated frequencies across 114 thermal tolerance alleles (B) per migration level are shown in transparent colors, and the mean is shown in bold.

Effects of genetic architecture on persistence predictions

Theoretical analyses suggest that the genetic architecture of adaptation will affect extinction probability during rapid environmental change (19). Our data allow us to ask whether the genetic architecture of climate adaptation affects predictions for climate change outcomes. We examined how the number of loci involved in thermal tolerance changed population trajectories under future warming scenarios. Extinction time was faster when more loci were needed to generate thermal tolerance (Fig. 4). For example, when there were 100 to 200 adaptive loci, time to extinction was 40 to 50 years under the RCP8.5 conditions. But with <50 adaptive loci, time to extinction was 60 to 80 years. Above 400 loci, simulations using RCP4.5 no longer led to persistence, and populations went extinct, showing that genetic architecture affects both the rate of decline and the ultimate outcome. These results are also reflected in the amount of time populations spent at low population sizes (below 50% of initial size) (Fig. 4). This simulation provides rough guidelines on the relative importance of genetic architecture: When thermal tolerance is controlled by fewer loci, adaptive change can be faster, and populations lose fewer individuals to natural selection.

Fig. 4 Impact of genetic architecture of thermal tolerance on population persistence under climate change scenarios.

Time to extinction (A and B) and total time (<50%) of initial population size (C and D) are shown. (A and C) Results for simulations conducted in increments of 10 between 10 and 200 loci. (B and D) Results show a broader sampling, in increments of 100 between 10 and 1000. Lines show means across 10 simulations per scenario.

Ecological determinants of persistence

Ecological parameters also affect evolutionary outcomes under warming scenarios. We tested the impacts of intrinsic population growth rate and initial population size on future persistence (Fig. 5). Growth rate (rmax) greatly influenced population dynamics: Populations with lower growth rates were more likely to go extinct and went extinct more rapidly. On the other hand, moderate changes in population size did not affect extinction estimates: Simulated populations from 500 to 3000 showed similar patterns across all warming scenarios.

Fig. 5 Heat map showing extinction times across a range of population growth rates and initial population sizes for all RCP scenarios.

Color reflects mean time to extinction for 10 replicate simulated populations. White cells indicate that none of the simulated populations experienced extinction.


Our simulations of one high-latitude coral population show that the rate and limits of adaptation to temperature may be sufficient to prevent the extinction of some populations under mild future emissions scenarios. On the basis of the parameters in our models, increased temperature leads to selection for higher frequencies of warm-adapted alleles. When increases in temperature are slow enough, the demographic deficit caused by selection is made up by increased population fitness and higher population growth rates. By contrast, in our model, more severe emissions outpace the evolutionary response, and mortality from selection is high enough to drive population extinction. The tipping point for the switch from adaptation to extinction as temperature change quickens is likely to be system-specific, but the basic scenario we see in our model is likely to be more general.

Our finding that faster warming in the future will likely decrease the probability of persistence is also seen in experimental studies of response to environmental stress in model systems such as Escherichia coli (20) and yeast (21). In these cases, rapid environmental change leads to adaptive collapse: Marked demographic decline as mortality increases, followed by increased loss of adaptive loci due to drift, and extinction before adaptation can restore higher fitness. In addition, models show that if the rate of environmental change outpaces adaptation, then populations become maladapted and unable to restore positive growth rates (20).

A striking feature of our results is that the cool-water population in Rarotonga harbors alleles identified as warm-adapted in populations located 800 km away. Although determining whether this is a more general phenomenon will require more research on more species, there are several reasons to think that there is wide occurrence of warm-adapted alleles in marine taxa. Many marine species have long-distance dispersing larvae, connecting populations over large distances (22, 23). These dispersal patterns tend to distribute alleles broadly, even as they allow buildup of frequency differences between locations. In other studies of adaptive divergence of long-distance dispersing marine populations, alleles are seldom private, and polymorphisms tend to be seen across large ranges (9, 24). As a result, clines across environmental gradients tend to be shallower in marine species, reflecting the balance between dispersal and selection (22, 23).

Major assumptions in the type of simulation we conduct here involve robust estimations of the allele frequencies of adaptive loci, the genetic architecture of thermal tolerance, and the map between environmental and genetic traits. A highly polygenic basis of thermal tolerance has been found in a number of systems (7, 8), including corals (15), suggesting that polygenic models are the most appropriate. However, some corals have specific loci with major effects on survival in the face of poor water quality (25) or sustained heat (26). For a precise estimation of genetic architecture, there is a need for further development of analytical tools, because current methods for detection of loci under selection can be riddled with false positives (27, 28). In addition, because genetic architecture estimated from natural populations can reflect environmental parameters unrelated to the key question, in our case temperature, further studies across replicate populations are needed for a more precise estimation of architecture. Careful consideration of experimental design and parameter choice is therefore necessary. We show that a rough estimation of genetic architecture might suffice for these types of predictions, because there is little difference in model outcome when more than ~100 loci encode the thermal tolerance trait. However, discovery that future climate tolerance was enabled by a few genes of major effect would substantially change the model.

Another area needing more information is the accurate estimation of phenotype, in our case, thermal tolerance, produced by different multilocus genotypes. We assume that there are many loci involved in heat tolerance and that all loci contribute equally without interactive effects, initial assumptions that are common in quantitative trait evolution. For example, in plant systems, numerous genetic variants predicted from genome scans across broad-scale contemporary environments can, together, predict performance under laboratory stress conditions (7, 29). These high-throughput measurements of tolerance and fitness are more difficult in nonmodel systems. Future efforts should therefore focus on precise and quantitative phenotypic measures to refine the map between genotype, phenotype, and future environment.

Future work should also focus on developing more realistic ecological scenarios within the simulation framework. Because we focused on a single isolated population, we used an extremely simple population genetic simulation with equal reproductive output from all surviving individuals, although fecundity can vary markedly and is influenced by factors such as depth and size (30, 31). Other extensions could include intermediate-latitude populations, integration of dispersal and connectivity, and even ecological factors such as resource availability and competition [reviewed by Hoban (32)] to create broader ecological frameworks that may be applicable to more species embedded in metapopulations.

Although our simulation provides a simple framework in which to begin using genomic data to predict climate change response, there are many, more complex interactions that would be beneficial to integrate into future models [reviewed by Bay et al. (33)]. The potential for plasticity to buffer short-term exposure to increased temperatures is well known for corals (14, 34, 35), but evidence in other systems shows that this plasticity can either promote or hinder adaptation (36). More careful work on dissecting the interaction between plasticity and adaptation of thermal tolerance limits in corals is needed to parameterize this model.

Another potential component in coral thermal tolerance is the role of the algal symbiont (genus Symbiodinium) (37, 38). Future studies could pair our model for host adaptation with local acquisition of adapted symbionts and the potential for symbiont shuffling in response to increased temperatures (39). Finally, alleles associated with increased thermal tolerance could be maladaptive in other conditions. For example, Howells et al. (13) found that warm-adapted individuals from the northern Great Barrier Reef experienced bleaching and mortality during winter low temperatures when transplanted onto the southern Great Barrier Reef. Future incorporation of these and other considerations will improve our ability to accurately predict population-level responses to ongoing climate change in other systems. Despite these caveats, simple adaptive models of coral populations in different future scenarios have the benefit of showing the potential value of genetic diversity in current populations and of the impact of relaxing assumptions that tolerance traits are unchangeable in climate envelope predictions.

Given the importance of adaptive landscapes in future climate responses in corals and other species, we used our integrated model to ask whether outplanting efforts could realistically increase the adaptive capacity of this high-latitude population. Adding a few corals with highly adaptive gene combinations could sometimes rescue the populations in our simulations. However, it is important to put the scale of outplanting in context. The level of assisted gene flow we simulated (10 individuals) is equivalent to replanting 1% of the population with adult, reproductive colonies every year. Accounting for expected levels of transplant mortality (17) and death of juveniles, this would require transplanting nearly 2 to 5% of the population annually until the year 2100. A bigger population than the small 1000 individual size we assume here would require proportionately more input. This would be a substantial annual effort with ongoing financial support for nearly a century.

Although these simulations suggest that this effort could have important adaptive consequences, it is important to note that incoming individuals must be genetically preadapted to future scenarios: Adding large numbers of individuals from populations with few heat tolerance alleles did not affect the outcome. This suggests that extreme care must be taken when selecting individuals for outplanting. Logistic and precautionary considerations (such as introduction of invasive species or disease) strongly suggest that only local outplanting should be conducted among local reefs with different thermal adaptations. Overall, there are many scientific, legal, and ethical concerns that should be carefully considered for any relocation project (40, 41).

A major outstanding question in climate change biology is whether the rate of adaptation will keep up with the rate of environmental change. We show that the current adaptive inventory of one coral population in Rarotonga is sufficient for moderate but not strong climate change. A growing number of studies show standing genetic variation in climate-related traits in many species, often involving large numbers of loci (7, 10, 24, 42). Other studies have become much better at predicting patterns of future environmental change on a local level. It may now be broadly possible to combine genomics and environmental modeling into a predictive framework that asks when and where evolution will be a major force in climate response.


Experimental design

Transcriptome sequencing. We collected fragments from 30 individual A. hyacinthus colonies from the back-reef on the southern shore of Rarotonga, Cook Islands. The fragments were part of a physiology experiment and, as such, had been subjected to experimental acclimation and heat stress or control conditions before preservation. Each fragment was then placed in RNA-stabilizing solution and stored for 24 hours at 4°C, after which it was stored for 1 to 2 weeks at −20°C. Samples were transported to Stanford University at room temperature and stored at −80°C until extraction. Total RNA was extracted from each fragment using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions, and complementary DNA libraries were prepared using the Illumina TruSeq RNA Library Prep Kit v2. Samples were randomized and, along with other samples from the same species, barcoded and multiplexed to sequence 12 samples per 50–base pair (bp) single-end lane. Pooled libraries were sequenced on an Illumina HiSeq 2000 at the Huntsman Cancer Institute at the University of Utah.

Statistical analysis

SNP identification. Adapter sequences and low-quality (Q < 20) regions were trimmed from raw 50-bp reads, after which short (<20 bp) reads were discarded. Quality-filtered reads were then mapped with bowtie2 (43) to the existing reference transcriptome for A. hyacinthus (44). After discarding duplicate reads, SNPs were identified using the Genome Analysis Toolkit (45), following the best practices published by the Broad Institute.

Because fragmentation in corals can result in multiple colonies of the same genotype or “genet” and skew population genetic analyses, we used hierarchical clustering to identify and collapse samples with very high genetic similarity (<0.1% divergence). SNPs that could be genotyped in all samples were used to build a distance matrix, and for pairs of samples with fewer than 10 sequence differences of this set of 9944 SNPs, we retained the sample with the largest number of genotyped SNPs. Downstream analysis was then conducted using only one representative from each genet.

Variation in thermal tolerance–associated SNPs. Previously, we identified SNPs associated with thermal tolerance across the A. hyacinthus transcriptome (11). For this, we took advantage of a thermal gradient on Ofu Island in the National Park of American Samoa that exists over a very small spatial scale. Two adjacent back-reef pools are distinct in their thermal variation: The HV pool occasionally exceeds 35°C, whereas less than 1 km away, the MV pool rarely exceeds 32°C. Acute heat-stress experiments confirm that A. hyacinthus colonies from the HV pool have higher thermal tolerance, even after a year of common garden (14). Using transcriptome sequencing, we genotyped 23 A. hyacinthus individuals from the two pools at 15,399 SNPs across the entire transcriptome. A combination of bootstrapped FST outlier and individual genotype by environment correlations identified 114 SNPs with strong evidence for selection associated with thermal environment (table S1).

Samples from the previous study (11) were combined with the new Rarotonga data set to examine patterns of genetic variation across the transcriptome. Quality filters were applied to the entire SNP set, discarding SNPs with fewer than six samples genotyped per population, as well as nonbiallelic SNPs and singletons. We calculated population differentiation between the two Ofu pools and between Ofu and Rarotonga using the Weir and Cockerham FST implemented in the hierfstat package in R. A principal components analysis was calculated for the entire set of SNPs and the subset of 114 SNPs associated with thermal tolerance using the prcomp function in R.

Simulations of future population size. On the basis of genetic parameters from the SNP data and demographic parameters derived from the literature, we used a simple logistic model for density dependence coupled with a genotype-specific survival function to simulate temporal changes in population size under future climate scenarios. The population size of the corals at year i + 1 was based on the number of colonies that died the previous year (Di), recruitment (Ri), and the standing stock of colonies (Ni)Embedded Image

The annual survival probability of each colony depend on the proportion of thermal tolerance alleles—those more common in the Ofu HV pool (11)—at 114 thermal tolerance SNP loci possessed by that colony compared to the optimum number of thermal tolerance alleles that were required for the thermal environment occurring that year. Initial diploid genotypes of each colony at all 114 thermal tolerance loci were simulated by binomial draws based on observed allele frequencies at each locus in Rarotonga. In subsequent generations, genotypes were simulated on the basis of allele frequencies at each locus of survivors from the previous generation. Survival of colony j at year i was modeled as a Gaussian function of the difference between the fraction of thermal alleles in colony j, Tprop(j), and the local optimum at year i, Topt(i). The SD (σ) of the Gaussian function was varied between 0.1 and 0.2Embedded Image

At each yearly interval, we determined the survival or loss of each colony by binomial draw based on the probability from this function. We then summed the number of colonies that perished as Di. Finally, we recorded the frequencies of thermal alleles at each locus for each of the NiDi colonies that remained in the population at year i.

In our simulations, all individuals had equal probability of reproduction. Recruitment depends on the number of surviving colonies, the intrinsic rate of increase (rmax), and the carrying capacity (K)Embedded Image

Population growth rate (rmax = 0.38) was derived from a study measuring recolonization after severe bleaching, from which we calculated the average annual population growth rate for Acroporid species over the decade after the event (18). Effectively, this rate would allow each colony to produce 0.38 daughter colonies each year in the absence of any density-dependent effects.

The SD of the Gaussian survival function determines how steeply survival declines for colonies without the optimal number of thermal tolerance alleles. We determined the SD by a preliminary set of simulations. Here, we assumed that the sampled population in Rarotonga was adapted to the current environment there. We therefore simulated a scenario without warming [Topt = Tprop(Rarotonga)] for 100 reproductive cycles based on the number of thermal tolerance alleles per colony observed in the Rarotonga population. The appropriate SD of ProbSurv(i,j) should maintain both population size and genetic variance in the absence of changing selection pressure. We simulated a range of SDs from 0.1 to 0.2 and found that an SD value of 0.16 best maintained both population size and genetic variance (fig. S4).

At the end of each year, we created the distribution of the number of thermal tolerance alleles based on survivors from the previous year (in the NiDi colonies). Therefore, the allele frequencies within recruits were assumed not to change this distribution of allele frequencies. These allele frequencies were then used to simulate diploid genotypes for the next generation.

Our simplistic model was a modification of the basic Wright-Fisher model commonly used in population genetic analysis. However, it does not account for inheritance. We also ran simulations with a fully individual-based model, where parental genotypes across thermal tolerance loci were retained between years and offspring were generated from random pairings of existing colonies. Because we know little about recombination rates between loci without a genome or linkage map, we assumed no linkage between thermal tolerance loci. The details of this model, in the standardized format provided by Grimm et al. (46), can be found in appendix S1.

Warming scenarios. We used the SST model output for the Rarotonga grid cell (−159.5, −20.88) from the NOAA GFDL ESM2M for the suite of RCPs developed for the IPCC AR5 (47, 48). ESM2M has moderate climate sensitivity among coupled ocean-atmosphere general circulation models developed for AR5 and therefore provides a moderate estimate of future warming (47). We also applied a correction to the SST output so that the mean and variance give peak temperatures and decadal variability statistics comparable with HadISST1 observations (4, 49), variables that would otherwise overestimate exposure to heat stress in corals. We used the temperature of the warmest month each year for our simulations. To calculate change in temperature and therefore selection pressure, we used a baseline equal to the average of the warmest months over the 10 years before sampling (2005 to 2015). We also included a “control” scenario with no warming. For each RCP scenario, we ran 10 replicate simulations. To determine whether using a different climate model would affect our results, we also ran a set of simulations on mean outputs across 35 IPCC AR5 models.

Impact of warming on optimum genotypes. To convert changes in SST to changes in Topt, we used genetic and environmental data for the MV and HV pools on Ofu based on 2 years of temperature data (2010 to 2012). HOBO loggers were placed at A. hyacinthus colonies and recorded the temperature every 10 min. We used temperature records from five individuals per pool, representing a total of 417,475 temperature readings per pool. We assumed that selection is the result of the highest temperatures in this system, because mean temperatures are approximately equal (16). We calculated the temperature extremes reached in each pool as the 1% highest temperature values recorded across all five individuals (T99). For the MV pool, this value was 30.36°C, whereas for the HV pool, it was 31.37°C. We then calculated a scaling factor (S) between degrees Celsius SST and Topt, which converted changes in temperature for each RCP scenario to changes in Topt at year i. We also assumed that changes in SST each year i (ΔSSTi) mirror changes in the 1% highest temperatures seen over timeEmbedded ImageEmbedded Image

After calculating S, we tested the sensitivity of the simulations to a range of scaling factors (fig. S5).

Sensitivity to genetic architecture and ecological parameters. The primary simulation was based on our best-estimated genetic and ecological parameters; that is, we used observed frequencies across 114 SNPs identified in our previous study (11), a realistic estimate of population size and a literature estimate of intrinsic growth rate. To explore the impact of genetic architecture of the thermal tolerance trait (in this case, the number of loci involved) on predicted persistence and extinction, we varied the number of loci simulated, randomly sampling from 114 SNPs with replacement. We conducted both a broad-scale examination, with the number of loci ranging from 10 to 1000, and a fine-scale examination, with the number of loci ranging from 10 to 200. These were simulated under all four RCP scenarios. Here, we simulated 10 replicates for each combination of parameters and calculated two metrics: (i) time to extinction and (ii) time below 50%—the mean number of years in which the population was less than half the initial size. Because corals are ecosystem builders and providers of coastal protection, the time below 50% metric provided a rough estimate of the persistence of these benefits.

We also examined the sensitivity of our simulations to ecological parameters: intrinsic population growth rate and initial population size. We varied intrinsic growth rate (rmax) from 0 to 1 and initial population size from 500 to 3000. For each value, we simulated 10 replicates and calculated the mean time to extinction across all replicates.

Impact of assisted gene flow. We conducted additional simulations to test whether migration of heat-resistant corals, for example, through transplantation from a more tolerant population, into the Rarotonga population could prevent extinction. For this, we simulated RCP6.0 with heat tolerance determined by 114 SNPs and a range of migration levels (0, 10, 25, or 50 individuals/year). Each year, we simulated individuals on the basis of allele frequencies from the HV pool and added them to the simulated population before selection was applied. We also conducted a control to rule out a strictly demographic effect where added individuals were simulated on the basis of the frequencies from the initial Rarotonga population. For each level of migration, we simulated 10 population trajectories and corresponding allele frequencies at each generation.


Supplementary material for this article is available at

appendix S1. Description of the individual-based model.

fig. S1. Comparison of probabilistic and individual-based models.

fig. S2. Sensitivity to choice of climate model.

fig. S3. Demographic impacts of migration.

fig. S4. Results of simulations for determining the width of the survival function.

fig. S5. Effects of the chosen quantile for the calculation of scaling factor (S) between degrees Celsius SST and Tprop.

table S1. Source data.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge members of the Palumbi Lab for assistance with field collection and useful discussion on analysis. Funding: This project was funded by the Stanford University Office of International Affairs; the Stanford Center for Computational, Evolutionary, and Human Genomics; and the Gordon and Betty Moore Foundation. Author contributions: R.A.B., N.H.R., and S.R.P. designed the experiment and collected the samples. R.A.B. conducted the laboratory work. All authors analyzed the data and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Genotypes for candidate SNPs are provided in table S1. Reads for RNA-seq data are uploaded to the National Center for Biotechnology Information Sequence Read Archive (BioProject ID: PRJNA411943).

Stay Connected to Science Advances

Navigate This Article