UC Frontiers of Biogeography Legacies of an ice-age world may explain the contemporary biogeographical provinces of corals

This study hindcast the geographic distribution of 18 Indo-Pacific scleractinian coral species, with different sensitivities to modern heat stress, into the last glacial maximum (LGM), some 18,000 years ago, when sea-surface temperatures were 2–4 o C cooler and sea level was ~130 m lower than on contemporary reefs. Identifying geographic provinces from the past may provide clues into genetic affiliations through time and provide some insight into how some coral species might respond to contemporary climate change. Coral habitat in the Indo-Pacific was reduced by 70% in the LGM. We identified five Indo-Pacific biogeographical provinces for corals during the LGM — (i) the western Indian Ocean, (ii) Southeast Asia, (iii) Indonesia and northwestern Australia, (iv) northeastern Australia and the Pacific Islands, and (v) the eastern Pacific. These provinces align with provinces recently identified using genetic markers. Given that the Quaternary was dominated by glacial conditions, the distributional legacies left through glacial dominance over that 2.6 million-year period may have also left evolutionary legacies. These glacial legacies may explain why contemporary corals live close to their upper thermal thresholds, which in turn have major consequences as we move into unprecedented rates of ocean warming. species distribution modeling, refugia, last glacial maximum, biogeographical provinces.


Introduction
For over 100 years biologists have contemplated the influence of glacial-interglacial cycles on terrestrial and marine populations (Agassiz andBettannier 1840, Reid 1899). Twenty-three glacial-interglacial cycles have occurred in the Quaternary (within the last 2.6 Ma) (Lisiecki and Raymo 2007), and the Earth has been in glacial conditions for 80% of the time over the last 700,000 years (Colinvaux 1987). During glaciations, tropical sea surface temperatures (SST) dropped by at least 2-4 o C (Guilderson et al. 1994, Koutavas et al. 2002, Lea et al. 2006) (supplementary Fig. S1) and sea level fell ~130 m (Chappell andShackleton 1986, Lambeck et al. 2014). These changes caused major reductions in shallow marine habitats and changes in coastal marine life (Wise and Schopf 1981, Valentine and Jablonski 1991, Lisiecki and Raymo 2005, Norris and Hull 2012. Recognizing that the Holocene period, which began approximately 12,000 years ago, is an anomalously warm time within the Quaternary, well above the average temperature of the last 2.6 Ma, it is important from an evolutionary perspective not to think of glacial conditions as differences from today, but also as differences from the preceding glacial period. It is also important to consider to what extent glacial periods have shaped distributional legacies and have contributed to biogeographic patterns (Hewitt 2000, Pellissier et al. 2014).
The trajectory of cooling, following the peak of Miocene warmth c. 18-14 Ma, was relatively gradual until the onset of Quaternary glaciations (Fig. S1). Beginning about 2.6 Ma a recognizable pattern of ~40,000 year glacial cycles ratcheted SSTs downward over the next 2 Ma (Lisiecki and Raymo 2007). Between 1 Ma and 0.5 Ma the glacial periods became stronger and longer, such that between 0.5 Ma and modern times there have been five glacial-interglacial events, each lasting ~0.1 Ma. The events were notably thermally asymmetric, with an 80,000 year cooling trend, reaching a glacial maximum, giving way to very rapid interglacial warming that typically lasted 5000-10,000 years. The downward path of temperatures was interspersed with precessional (19,000-23,000 year periodicity) modulations that provided a pause in cooling rather than true thermal reversals in tropical ocean systems. Thus, the physiologically novel conditions for marine life during the glacial-interglacial cycles would have followed the glacial maxima when abrupt warming during deglaciation may have acted as a significant stressor. Atmospheric and oceanic circulation patterns were also weaker during glacial periods than currently (Rahmstorf 2002), resulting in reduced genetic connectivity during the glaciations (Norris and Hull 2012). Such reduced connectivity may have led to reduced genetic diversity (Hewitt 2000). However, these events may have also spurred local adaptation and speciation events (Wilson et al. 2009), leaving distributional legacies that extended into interglacial periods (Potts 1983, Benzie 1999, Hewitt 2000.
Studies have identified distributional legacies by identifying biogeographic provinces in which marine fauna are closely related. Although early work relied heavily on taxonomic affinities (Veron 1995, Spalding et al. 2007, Kulbicki et al. 2013, early genetic analyses of giant clams  and Linkia sea stars  showed biogeographic provinces and a clear decoupling between genetic affinities and contemporary oceanographic currents, particularly in the western Pacific Ocean. Benzie (1999) showed that the gene flow trend of three species of giant clams ran north-west to south-east, which was perpendicular to the contemporary South Equatorial Current that runs north-east to south-west through the region (Benzie 1999). This mismatch shows that past sea levels and past episodes of dispersal had lasting effects on modern marine species (Ludt and Rocha 2015). Recently, molecular analyses of 56 marine taxa revealed five major biogeographic provinces across the Indo-Pacific (Crandall et al. 2019) (supplementary Fig. S2). However, it is still unclear whether the entire coral triangle, situated at the intersection of the Indian and Pacific Oceans, is one biogeographic province or is comprised of many subregions with similar genetic affinities (Barber et al. 2000, Briggs and Bowen 2012, Keith et al. 2013. To gain insight into the major coral provinces of the LGM, we hindcast the geographic distributions of representative scleractinian coral species. Scleractinian corals first appeared in the Triassic, some 237 million years ago, after the Permian extinction (Stanley and Fautin 2001).
However, modern reef corals did not really diversify until the mid-Jurassic period, with the simultaneous emergence of their symbiont dinoflagellates (family Symbiodiniaceae) (LaJeunesse 2020). Corals faced numerous environmental challenges through the Cretaceous-Tertiary extinction event and through the Eocene thermal maxima (Veron 1995, Budd andJohnson 1999). Reef-building corals also face serious challenges today because of thermal-stress events (Hughes et al. 2018). Indeed, scleractinian corals have become sentinel species of modern coral reefs because of their sensitivity to ocean warming. Climate change is causing marine heatwaves that induce corals to bleach (Sully et al. 2019). Coral bleaching is a result of the breakdown of the symbiosis between scleractinian corals and zooxanthellae-their photosynthetic symbionts. This study hindcasts the geographic distributions of 18 Indo-Pacific scleractinian coral species, with different sensitivities to modern heat stress, into the LGM using past SST and shelf-exposure conditions. Identifying geographic provinces of scleractinian corals from the past may provide both clues into processes of genetic differentiation through time and insight into how some coral species might respond to contemporary climate change.

Biological data
Species distribution models were implemented to compare contemporary distributions of Indo-Pacific corals with hindcast predictions for the LGM. Species distribution models associate the geographic distribution of a species with broad-scale environmental conditions, combining the geographical space in which a species of interest is found with its theoretical environmental space. The environmental space is used to predict the niche of the species, and then that environmental space is projected back onto geographic space, for which environmental conditions are known (Elith and Leathwick 2009, Cacciapaglia and van Woesik 2015, Cacciapaglia and van Woesik 2018. Species distribution models were run for contemporary corals and for hindcast corals using environmental conditions in the LGM. The models predicted the distribution of 18 Indo-Pacific coral species, selected based on their short-term response to heat stress (Loya et al. 2001). Six of the selected corals are considered heat-stress tolerant species: 1) Goniastrea aspera, 2) Porites lobata, 3) Cyphastrea chalcidicum, 4) Dipsastraea speciosa, 5) Favites halicora, and 6) Leptastrea pruinosa. Six of the selected coral species are considered heat-stress intolerant species: 7) Porites horizontalata, 8) Montipora aequituberculata, 9) Acropora hyacinthus, 10) Acropora digitifera, 11) Seriatopora hystrix, and 12) Stylophora pistillata. Six coral species were considered to have intermediate heat-stress tolerance, 13) Acropora cerealis, 14) Dipsastraea pallida, 15) Favites pentagona, 16) Galaxea fascicularis, 17) Goniastrea retiformis, and 18) Platygyra sinensis.

Environmental parameters
SST range and SST minimum for contemporary climate were derived from Aqua-MODIS (Moderate Resolution Imaging Spectroradiometer) at a 9.2 km 2 resolution from 2002-2009 to explore niche limits. These data were collected from BioORACLE (Ocean Rasters for Analysis of Climate and Environment) (http://www.oracle.ugent.be/) (Tyberghein et al. 2012). SST data for the LGM were collected from GISS CMIP5 projection as the range and minimum, at a resolution of 2.5°. These data were downscaled to fit the resolution of contemporary SSTs.
Contemporary photosynthetically active radiation (PAR) data were derived from SeaWiFS (Seaviewing Wide Field-of-view Sensor) observations at 9.2 km 2 resolution from 1998-2007. A global raster of PAR during the LGM was reconstructed using the 'palinsol' package (Crucifix 2016) in R (R Core Team 2018) using current PAR as a template ( Figure S3). Turbidity was derived from SeaWiFS data using (Kd490) and PAR to calculate the irradiance at a 3 m depth (Cacciapaglia and van Woesik 2016) with the following equations derived from Gallegos and Moore (2000) and Pierson et al. (2008): where KdPAR is the diffuse coefficient of PAR and Kd490 is the diffuse attenuation coefficient of downwelling spectral irradiance wavelength 490 nm. Irradzm is the proportion of irradiance reaching the depth Zmax. This equation was applied to mean PAR both during the contemporary timeframe and the LGM. Turbidity data were binned to prevent no-analog estimates of a random intercept, using Akaike's Information Criteria to find the most parsimonious bin size. Turbidity was resampled for the LGM and it was assumed that tidal amplitude would follow similar global patterns in the past. Changes in sea-level and the shape of the continental shelf during the LGM would have caused local patterns of turbidity to change -this was partially mitigated by resampling to a courser 5 o resolution to match the LGM datasets so that geographic averages could be observed at a broad scale, rather than using fine scale data that would be inaccurate along a shifting coastline.

Mechanistic constraints
Five mechanistic constraints, or 'masks', were imposed on coral distributions based on physiological thresholds to increase model parsimony and impose realistic constraints. First, a 6 bathymetric depth mask was imposed to exclude habitat ≥ 30 m. This mask changed between the two different time periods -to include 0-30 m for contemporary coral reefs and 130-160 m during the LGM. This bathymetric mask was obtained using NOAA bathymetric data from the getNOAAbathy function in the R package 'marmap' (Pante and Simon-Bouhet 2013). This first mask neglects tectonic rebound from glacial weighting as well as potential reef growth. Second, a river outflow mask was imposed to prevent false positives in freshwater outflow of 12 major tropical rivers. The size of each river outflow mask was based on the drainage size of each major river (Cacciapaglia and van Woesik 2015). The position of the river delta was adjusted during the LGM projection by identifying the geographic carvings on the continental shelf, using Google Earth, when sea level was around 130 m lower than present. We did not however consider changes in drainage areas and precipitation with lowered sea level. Third, a dispersal mask was used to mask out the Atlantic Ocean, because the Isthmus of Panama has acted as a physical barrier to dispersal in both the LGM and currently. Fourth, a thermal minimum mask was imposed to remove habitat of ≤18 o C, to prevent false positives in locations that were too cold to support reef-building corals both on contemporary reefs and during the LGM. The LGM <18 o C mask was also used to determine how much available habitat was precluded by this thermal minimum during the LGM. Fifth, a limitation in minimum irradiance mask was imposed following Kleypas et al. (1999), where coral presence was masked anywhere ≤ 21.6 μEm -2 d -1 .
The geographic extent of the study was limited to between 37 o N and 37 o S.

Model implementation
A Bayesian generalized-linear-mixed-effects model from the R package 'blme' (Chung et al. 2013) was used to examine the response of coral species to climatic parameters following: , where yi is the observed binary response of a species (presence or absence) at location i, pi is the probability of a species being present at location i. The binomial distribution is the sum of the Bernoulli distributed cases of presence or absence of species considered for all locations n. Presence points were sampled from the known species distributions using Veron (2000), IUCN's Red List http://www.iucnredlist.org/technical-documents/spatial-data, and Wallace (1999) on reef locations derived from ReefBase. Absence points were sampled using a probability of detection algorithm following Cacciapaglia and van Woesik (2015) where absence points were sampled from reefs in ecoregions that had >95% probability that the given species was absent. The number of presence and absence points that were sampled followed a ratio depending on the 7 reef area within occupied and unoccupied ecoregions. The coefficient βo is the intercept of the model, β1, β2, and β3 are slopes that define the extent of the relationship between the response variable and the environmental covariates, xi (range of SST) and zi (range of PAR) at site i, and here turbidity, a, was expressed as a random-error effect, which followed a gamma distribution. We used the range in SST and PAR to inform our model of niche space. We considered that range was a better informant of niche space than averages, and range was more parsimonious than using both minimum and maximum. Fixed effects were allocated to SST range and PAR range, and a random effect was allocated to turbidity. This model was trained using contemporary environmental parameters, and hindcast for environmental parameters reconstructed from the LGM. A comparison of the model outputs was used to determine the changes in coral species distributions between current and those hindcast for the LGM. We tested various other models, although the Bayesian Generalized Linear mixed models outperformed Generalized Linear Models (GLM), Generalized Additive Models (GAM), and bioclim models using the AUC (Area Under the receiving operating Curve) scores. The only model to outperform the GLM was a RandomForest that scored a higher AUC in a pixelated hitor-miss fashion, which did not represent a biological system where regional patterns are more realistic. These models were all evaluated in Cacciapaglia and van Woesik (2015).

Model evaluation
A k-fold procedure in R (R Core Team 2018) was used to test and validate the model output. The k-folding used 80% of data points for training and the remaining 20% of the data points for testing. The Area Under the 'Receiving Operating Characteristic' Curve (AUC) was used to evaluate model integrity, where ≥ 0.8 was considered a robust model. This evaluation was repeated 100 times in each of the 25 model runs for each of the 18 coral species, ensuring that the highest scoring AUC models were used for the predictions and allowing for confidence in each projected species. The greatest sum of sensitivity and specificity was used to define the occupancy thresholds, which were converted to presence and absence maps.
A distance matrix was run on the coordinates gained from the projected distributions of the 18 scleractinian coral species for the conditions that were expected in the LGM. Reef habitats that were contiguous and ≥ 0.25° in geographic area when aggregated were used to calculate the distribution of major coral provinces. The distance matrix was run using the LGM aggregates in a pairwise comparison, finding the geographic distance between the locations of each major coral province to show the distance of contemporary coral reefs from major coral provinces that were hindcast for the LGM. The combined distribution for all 18 coral species was used to determine the aggregate when ≥ 12 of the 25 runs showed that a coral species was present in the model output. These aggregated reefs and their distance from one another were used to identify macrorefugia during the LGM. A hierarchical clustering analysis using the function 'hclust' in the R 'stats' package was run on the modern and LGM time periods to: 1) determine whether our results for modern corals match contemporary biogeographic regions, and 2) identify likely biogeographic regions during the LGM. A one-way analysis of variance (ANOVA) was run on the change in each species habitat occupancy between modern projections and LGM projections using the three coral-species groups as factors that were categorized as, (1) heat-stress intolerant,

Results
The models predicted that the 2-4 o C lower SST during the LGM shifted corals closer to the Equator and caused an estimated 0.75 % loss of coral habitat (~17,000 km 2 ). At the same time, the available shallow-water habitat for corals in the Indo-Pacific was reduced three-fold ( Figures  1 and 2, Figure S4, Table 1). This change in habitat availability was largely a function of the lack of available shallow habitat beyond the continental shelves. The change in minimum SST, which was responsible for a reduction in coral cover by an average of 17,000 km 2 among all species (table S1), had less impact than habitat change caused by the loss of shallow reef area during the LGM-a difference of 2.4 million km 2 of habitat above 30 m depth throughout the Indo-Pacific tropics.  LGM (see also the kmz files that project each species in the LGM and modern on Google Earth in Appendix S1/Earth)). The most diverse epicenters of coral occupancy were located in Southeast Asia, northeastern Australia, northwestern Australia, and the Pacific Islands during the LGM, with two secondary localities in the Seychelles and in the eastern Pacific (Figures 1, 2, and 3). Major differences between modern and LGM coral distributions were apparent in areas where the shallow-water habitat was available during the LGM, such as in the Maldives, Seychelles and Mauritius Islands (Figure 1; Figure S5). The smallest change in predicted distributions occurred in areas with steep continental shelves, as on most of the central coast of eastern Africa. In the eastern Pacific Ocean, there was more coral prevalence during the LGM, which diminished under modern climate and bathymetry (Figure 3). The French Polynesia region supported corals in the LGM ( Figure 3); see also the online kmz file that projects this figure on Google Earth), but did not provide a large shallow water area (> 0.25°), and therefore was not highlighted as an occupancy hotspot (Figure 3).  There was a positive relationship between the distributions of 18 contemporary Indo Pacific coral species and those hindcast for the LGMs ( Figure S6). The most prevalent contemporary Indo-Pacific coral species were also the most prevalent corals during the LGM. Many coral species were ubiquitous during the LGM, although some species were more widely spread than others. For example, Porites horizontalata lost relatively less habitat (40%) than other species during the LGM (Table 1). By contrast, Porites lobata lost 66% of its habitat during the LGM (Table 1). Similarly, Montipora aequituberculata and Stylophora pistillata lost 65% and 64% of their habitat during the LGM and were predicted to have had lower habitat occupancy during the LGM than most other coral species. The species that lost the most habitat was Favites pentagona, with 70% loss of habitat. An analysis of variance (ANOVA) run on the three coralspecies groups found no significant difference (p= 0.37) in the change in habitat occupancy between modern projections and LGM projections among the groups. Species with restricted geographical distributions had less informative models, and some models had an AUC score below 0.8, our threshold for a robust model. These models are indicated with an asterisk in Table  1. We identified five Indo-Pacific biogeographical provinces for corals during the LGM -(i) the western Indian Ocean, (ii) Southeast Asia, (iii) Indonesia and northwestern Australia, (iv) northeastern Australia and the Pacific Islands, and (v) the eastern Pacific (Figure 4; Figure S7).

Discussion
The last glacial maximum During the LGM, when sea level was ~ 130 m lower than today, vast shallow-water systems that support diverse coral were dry land. For example, the Sunda and Sahul shelves in Southeast Asia, which today are the epicenters of marine diversity, would have supported lowland forest. As sea level fell, most reefs withdrew to the edges of the continental shelves where a steep seafloor maintained a shallow-water habitat. The loss of coral habitat by lowered sea level was substantially greater than the loss caused by oceanic cooling (see supplementary document, Table S1), but together these factors induced a 70% loss of coral habitat during the LGM. Di Nezio et al. (2016) developed models to hindcast the oceanic processes during the LGM for the western Pacific warm pool. They found that, like today, simulated oceanographic currents had a high westward surface velocity and a strong countercurrent just north of the Equator, generating connections between the Pacific and Indian Oceans.
Our comparison between the distributions of 18 contemporary Indo-Pacific scleractinian coral species and those hindcast for the LGM shows that, despite the loss of habitat during the LGM, the corals had a similar distribution to modern corals, particularly those corals in Southeast Asia and northeastern and northwestern Australia (Figure 2). With the exception of the heat-stress intolerant corals Stylophora pistillata, Montipora aequituberculata, and Dipsastraea speciosa, and some species with intermediate heat-stress tolerance, such as Platygyra sinensis and Favites pentagona, most corals could have occupied almost all available habitat inside the 18 o C constrained boundaries, indicating favorable conditions for equatorial corals during the LGM.
These model estimates agree with Kleypas (1997), whose analysis suggested favorable reef growth during the LGM near the Equator. This similarity was further supported by the ANOVA, run on the three groups in the present study (which were categorized based on heat-stress tolerance), which showed no significant difference in loss of habitat among the three coral groups during the LGM.

Biogeographic provinces
The present model identified 5 biogeographic provinces for corals during the last glacial maximum: (i) the western Indian Ocean, (ii) Southeast Asia, (iii) Indonesia and northwestern Australia, (iv) northeastern Australia and the Pacific Islands, and (v) the eastern Pacific ( Figure   4). These provinces would have formed the sources for deglacial expansion. At the onset of deglaciation, the rapid warming and sea-level rise allowed corals to radiate from biogeographic provinces to occupy newly available habitat on continental shelves. The two regions in Australia were likely further apart than the cluster analysis suggested (Figure 4), since during the LGM there was a landmass between the two regions and the temperature constraints do not allow for connectivity south of the continent. The northeastern Australia and the Pacific Islands province suggests high connectivity in the past on either side of the north-west to south-east, Australia-New Guinea axis. This province agrees well with genetic affinities (Benzie and and is decoupled from the modern South Equatorial Current, again suggesting past episodes of dispersal had lasting effects on modern marine species (Ludt andRocha 2015, Crandall et al. 2019).
Past literature showed that terrestrial migrational velocities would have had to be unrealistically high to keep pace with warming through the deglacial period, exceeding empirical observations by a factor of 10-100, which became known as Reid's Paradox . A solution to the paradox came from two different explanations. The first explanation argued for unusually long-distance dispersal events, amounting to jump-dispersal (Clark 1998). The second explanation argued for survival in microrefugia, far from the main populations (McLachlan and Clark 2005, Feurdean et al. 2013. In either case, founder effects and selection in small subpopulations may have shaped the genetics of isolated populations. When deglaciation conditions favored population expansions, these isolated populations, with distinct genotypes, formed local vanguards of a larger-scale invasion (Mosblech et al. 2011, Cheddadi et al. 2017).
Small isolated coral populations found in microrefugia have lower diversity, higher genetic load, and lower fitness than populations in macrorefugia (Mosblech et al. 2011). By extension, corals on reefs in Hawaii should be less genetically diverse than others, as they have been isolated from biogeographic provinces through glacial events. By contrast, coral populations on the Indonesian reefs are more likely to be genetically diverse because of their habitat stability through glacialinterglacial cycles. Indeed, Baums et al. (2012) showed that, for example, the heat-stress tolerant coral Porites lobata colonies located in the central Indo-Pacific had much higher allelic richness than the colonies found in the Hawaiian Islands. While these genetic results could be simply a consequence of geographical isolation (Selkoe and Toonen 2011), they could also reflect the location of the glacial provinces identified in the present study ( Figure 2).
The LGM provinces identified by our models for corals ( Figure 4) align with biogeographic provinces recently identified using molecular markers of 56 marine taxa (Crandall et al. 2019) ( Figure S2). Considering Crandall et al.'s (2019) results from west to east, they reported an Indian Ocean delineation east of the Seychelles, which aligns with our model results (Figure 4). They also showed a province that extends from western Indonesia into the central Pacific, near Samoa ( Figure S2). We predicted a similar province, but our models projected a delineation between Southeast Asia and Indonesia-northwestern Australia. Our delineation instead agrees with Barber et al.'s (2000) 'Marine Wallace's line' through the Coral Triangle ( Figure S2), which was predicted using molecular markers of mantis shrimp. Crandall et al. (2019) also found delineation west of the Hawaiian Islands and an eastern Indo-Pacific province that extended to the Marquesas Islands, although the latter province was not significantly different from their randomized models. These two delineations agree with our central Pacific delineation of the separate biogeographic provinces between the central and eastern Pacific Ocean (Figure 4).

Glacial legacies
Glacial periods have dominated glacial-interglacial cycles by a ratio of 4 to 1 (Colinvaux 1987) over the last 700,000 years, and therefore corals would have spent the majority of their recent evolution in glacial-type conditions. We hypothesize that these dominant glacial periods left evolutionary legacies and may explain why contemporary corals live close to their upper thermal thresholds, as cool temperatures were the norm through the Quaternary. In support, Schoepf et al. (2019) showed that Acropora aspera maintained heat tolerance under cooler temperatures but did not acclimate to warmer temperatures. Similarly, Sinclair et al. (2016) showed that typical ectotherms displayed a left-skewed (negative skewed) thermal performance curve, suggesting that modern ectotherms can adjust to cooler temperatures, but quickly reach maximum thermal limits. Pörtner (2012) also argued for an asymmetrical, or non-Gaussian, thermal niche, where optimal metabolic performance rapidly declines at high temperatures. We suspect that there is no reason that glacial traits were purged from populations during interglacial periods, and we should not be surprised therefore that Acropora is one of the most vulnerable species to ocean warming (Loya et al. 2001) when its phylogenetic expansion occurred primarily during glacial periods.
Our modern proximity to the upper thermal boundary of Quaternary history has been identified as a potential stressor (Colwell et al. 2008, Hoegh-Guldberg 2007, IPCC 2013. The most recent time that would have promoted selection for warmer-than-modern traits was at the peak of the last interglacial between c. 128,000 and 125,000 years ago, when SSTs were 1-1.5 o C warmer than contemporary SSTs (McCulloch and Esat 2000). The traits conferring upper levels of thermal tolerance were unlikely to have been maintained through the intervening ≥100,000 years when selection would have favored cool-adapted traits. Thus, the asymmetry of glacial temperatures dominating the Quaternary may explain why many corals, particularly acroporids, are sensitive to marine heatwaves, because scleractinian corals have not experienced such high SSTs for millennia. Persistent cool SSTs through the Quaternary may also help explain the extent of contemporary coral bleaching from thermal-stress events -which is widely considered the largest threat to contemporary loss of reef corals (Hughes et al. 2018).

Model refinement
Our work suggests that glacial biogeographic provinces may be just one factor in the varied processes that drive spatial differences in genetic affinities of corals. Nearly all climate reconstruction models, transfer functions, and predictive models necessitate the assumption that species remain in their modern-day niche through time. Some parts of the modeled responses to LGM condition are absolute (i.e., the requirement to be below sea level), whereas the thermal response could be more plastic than assumed. Recent studies on bleaching may shed light on the timeframe in which corals have the capacity to adjust. Sully et al. (2019) showed that Indo-Pacific corals may have adjusted to a 0.5 o C warming in SST over the last decade. This adjustment entails the complexities of acclimation, adaptation, and community shifts toward more temperature-tolerant coral species, although it does not demonstrate the capacity to escape an evolutionary thermal maximum. Additionally, Kenkel and Matz (2016) showed major differences in thermal tolerance within coral species that were no more than 10 km apart in the Florida Keys. These small-scale temporal and spatial differences suggest a greater need for highresolution spatial and temporal information on thermal tolerance of corals.
In conclusion, we suggest that the glacial dominance over the last 2.6 Ma have left distributional and evolutionary legacies of corals, which have consequences for contemporary corals close to their thermal thresholds, especially as we move into unprecedented rates of ocean warming. Indeed, contemporary corals are facing increased SSTs, which cause bleaching and mortality. Stenothermal species, with a narrow distribution and niche, will most likely fare poorly under future climate change, particularly if most of their distributions are outside of the area that supports high genetic diversity, such as outside of the Coral Triangle. Using a simple distance proxy from macrorefugia during the last glacial maximum, our study suggests that the highest diversity region to resist climate change is most likely located in the central Indo-Pacific rather than in the equatorial central Pacific.