Palau’s warmest reefs harbor a thermally tolerant coral lineage that thrives across different habitats

Ocean warming is killing corals, but heat-tolerant populations exist; if protected, they could replenish affected reefs naturally or through restoration. Palau’s Rock Islands experience chronically higher temperatures and extreme heatwaves, yet their diverse coral communities bleach less than those on Palau’s cooler outer reefs. Here, we combined genetic analyses, bleaching histories and growth rates of Porites cf. lobata colonies to identify thermally tolerant genotypes, map their distribution, and investigate potential growth trade-offs. We identied four P cf. lobata genetic lineages. On Palau’s outer reefs, a thermally sensitive lineage dominates. The Rock Islands harbor two lineages with enhanced thermal tolerance and no consistent growth trade-off. One of these lineages also occurs on several outer reefs. This suggests that the Rock Islands provide naturally tolerant larvae to neighboring areas. Finding and protecting such sources of thermally-tolerant corals is key to reef survival under 21st century climate change.


Introduction
Ocean warming with increased marine heatwave intensity is considered the most signi cant threat to coral reefs globally 1 . Seawater temperatures >1°C above historical summertime values can disrupt the coral-algae endosymbiotic relationship, leaving corals in a nutritionally compromised, pale state ("bleached") and vulnerable to death 1 . Since the early 1980's, millions of corals have died and thousands of acres of coral reef area have been lost as bleaching events have become more frequent, severe 2 , and widespread 3 .
As anthropogenic CO 2 emissions continue to rise, fueling further warming and more intense heatwaves, science and conservation efforts are increasingly focused on identifying thermally-tolerant coral communities that could survive ocean warming and potentially reseed impacted reefs either naturally or via restoration 4 . Environments with high variability are promising coral source sites because they may promote colony plasticity and harbor resilient communities [5][6][7] . In addition, habitats that currently experience conditions akin to those projected under climate change (e.g. higher temperatures, lower pH), are potential reservoirs of environmentally tolerant coral populations that may facilitate coral survival through the in ux of climate-adapted offspring (i.e. evolutionary rescue 8 ).
Despite warmer temperatures, there are more coral genera and higher coral cover (~60%) in the Rock Islands than on Palau's outer reefs 15 . Warmer temperatures and high coral diversity can facilitate selection for communities with higher thermal tolerance. During heatwaves associated with the El Niño Southern Oscillation (ENSO), the Rock Islands experience larger temperature anomalies than outer study is the rst to combine individual colony genetics with growth parameters and bleaching histories of responses to natural ENSO heatwaves through coupled genetic sampling and coral coring.

Results
Rock Island sites have consistently warmer temperatures and a higher diurnal range.
Multiple years of in situ temperature data from four Rock Island sites reveal these locations are chronically warmer than all outer reef sites (Fig. 1A). Across Palau's outer reefs, average water temperatures were 29.11°C +/-0.7, with a mean diurnal range of 0.32°C. Mean temperatures on Rock Island reefs were ~1.5°C warmer at 30.29°C +/-0.62, and with an 85% higher mean diurnal range of 0.59°C. Some days at Ngerchelong and Helen reefs (both outer reefs) had wide temperature ranges, but these were rare events compared to the more frequent high-range days at most other Rock Island reefs ( Fig. 1B). Taking the 90th percentile warmest temperatures each year, the top 10% of peak temperatures in the Rock Islands were above 31.08°C on average, compared with only 29.78°C on the outer reefs.
Porites cf. lobata in Palau form four distinct genetic lineages across habitats with varying temperatures, and lineages host the same symbiont types.
Microsatellite genotyping and RAD-sequencing analyses revealed that P. cf. lobata populations across the Palauan archipelago form four distinct lineages with high genetic divergence (Figs. 2-3; S1). STRUCTURE assignments of samples to genetic lineages were consistent between microsatellites and RAD-seq datasets, with 96% of samples assigning to the same lineage (Fig. S2). PCA and DAPC analyses supported the presence of four strongly differentiated lineages in RAD-seq ( Fig. 3B-C) and microsatellite data ( Fig. S1B-C). Please note, the color designation of lineages in gures is simply for visualization purposes and does not re ect differences in colony color.
As expected, F ST values between lineages were higher when calculated based on RAD-seq rather than multi-allelic microsatellite data, with ranges between 0.24 and 0.67 for RAD-seq data and 0.07 to 0.17 for microsatellite data (Figs. 3D; S1D). The ordering of pairwise differences remained the same between the two datasets. The dark blue (DB) and light blue (LB) lineages were the least differentiated (RAD F ST =0.24; microsatellite F ST =0.07), and had more admixed individuals. The high (>0.3) F ST values and the high lineage assignment probabilities suggested that some of these lineages could represent cryptic species of Porites lobata or other massive Porites species across the Palauan archipelago. The presence of several admixed individuals (Figs. 3A; S1A), however, indicated that these lineages likely still hybridize; though these signatures may also suggest past introgression. Here, we refer to these genetic groups as differentiated lineages, absent solid evidence of reproductive isolation or incompatibility that would support reclassi cation based on the biological species concept.
The four lineages were differentially distributed across Palau's reef habitats, with some lineages predominantly found in the warmer Rock Islands, while others were dominant on cooler outer reefs (Fig. 2). The DB lineage was predominant among outer reefs, while the LB lineage was the most widespread, though it occurred in higher proportion within warmer Rock Island sites (Fig. 2). The pink (PI) lineage was more common on outer reefs but represented a small fraction of the community at most sites, except in Ngerchelong (Fig. 2). The red (RD) lineage was entirely con ned to the Rock Islands, except one individual from Drop Off (Figs. 2; 3A). Within the Rock Islands, RD corals were common in Mecherchar and Risong and rare elsewhere (Fig. 2).
We found no differences in the symbiont composition of colonies using denaturing gel electrophoresis (DGGE) and ITS2 sequencing (Table S1). In fact, all colonies harbored Cladocopium (C15) symbionts, which are strongly associated with massive and branching Porites corals across the Info-Paci c, including Palau 16, 25,26 .
Lineages had different growth rates and thermal tolerances.
Skeletal density, calci cation, and linear extension rates differed among lineages ( Fig. 4A-C). The Rock Island-associated RD lineage had lower skeletal density than all other lineages (ANOVA, F=9.04, df=3; post-hoc Tukey test, p<0.05; Fig. 4A), as well as lower calci cation and extensions rates than the DB and LB lineages (ANOVA, F=5.96, F=3.10, respectively, df=3; post-hoc Tukey tests, p<0.05; Fig. 4B-C). The outer reef-associated PI lineage did not have signi cantly lower growth rates for any metric compared with the DB and LB lineages, but it was also not signi cantly different from the RD linage in calci cation or extension rates, showing an intermediate phenotype between the DB and LB lineages and RD ( Fig. 4A-C).
The presence or absence of high-density stress bands in mounding Porites cores provides valuable insight to differences in thermal tolerance among individuals and across habitats (e.g. Fig. S3). The widespread outer reef DB lineage showed the highest prevalence of stress bands (68%) during the 1998 ENSO, suggesting this lineage has low thermal tolerance. In contrast, the Rock Island-associated LB and RD lineages showed signi cantly lower stress band prevalence at 22% and 25%, respectively (X 2 =12.35, df=3, p-value=0.006; Fig. 4D). The PI lineage also showed low stress band prevalence, though we had the fewest number of samples from this lineage (N=6). During the less severe 2010 ENSO, lineages showed overall lower prevalence of stress bands across all lineages (Fig. S4). Notably, LB and DB had similar growth rates and skeletal density despite their differences in thermal tolerance in 1998 (Fig. 4). The RD lineage demonstrated high thermal tolerance and low skeletal growth and density, but these corals are con ned two of the lowest pH Rock Island sites 14,15 , which can slow growth and reduce density 27,28 (Figs. 2;4).

Discussion
Here we have combined, for the rst time, historical growth and bleaching responses to ENSO heatwaves in the eld with genetic data, to nd that Palau's Porites cf. lobata populations form distinct genetic lineages with differing thermal tolerances (Figs. 2-4). The unique environmental conditions of Palau's Rock Islands have likely promoted the development of thermally tolerant lineages (Fig. 1). Lineages most common in the Rock Islands can tolerate higher temperatures before bleaching and bleach less frequently (Fig. 4). When these lineages are found outside of Rock Island habitats, they still maintain higher tolerance, suggesting this is at least partially driven by genetic factors (Fig. 5D). Furthermore, the Rock Island-associated LB lineage maintains higher thermal tolerance across habitats, without showing any corresponding trade-offs in growth metrics, indicating this lineage may be particularly well suited for restoration efforts (Fig. 5A-C). Our results expand prior work examining the in uence of warmer and/or more variable habitats in promoting thermal tolerance (e.g. the back reef pools of Samoa 6 ; the Florida reef tract 7,29 , and Australia 30,31 ) to show that these patterns are repeatable across multiple reefs. Our results also reinforce the notion that diverse habitats can yield climate change resistant corals 32 .
Temperatures in the Rock Islands are both warmer and more variable (Fig. 1), which could enhance the thermal tolerance of resident corals in several ways. First, the higher average temperatures can serve as a lter, selecting thermally tolerant corals through their early life stages, including selecting against less thermally tolerant larvae that may arrive from outer reefs. Limited connectivity between the Rock Island and outer reefs as suggested by hydrodynamic model output [33][34][35] likely facilitates local adaptation, which is corroborated by the population differentiation seen here (Figs. 2-3). We saw strong genetic differentiation between lineages and found that reefs within the Rock Islands have different lineage compositions than outer reefs (Figs. 2-3; S1). We further showed that lineages differ in their growth rates and thermal tolerance (Fig. 4). The outer reef-associated DB lineage showed signi cantly higher bleaching (lower thermal tolerance) than other lineages during the 1998 ENSO event, while the LB and RD lineages, most common in the Rock Islands, had signi cantly lower bleaching (Fig. 4D). These patterns are consistent with the notion that Rock Island environments act as selective lters to produce lineages that are adapted to tolerate warmer temperatures. This is further supported by the LB lineage maintaining low bleaching prevalence in 1998 when living on the outer reefs, which points to a genetic basis for its thermal tolerance (Fig. 5D).
Second, daily temperature ranges between 0.5°C to 5°C have been found to be a better predictor of bleaching patterns than other environmental parameters 5 , and in both natural and experimental studies corals from variable environments often show increased thermal tolerance 6,31,36 . For instance, in American Samoa, Acropora hyacinthus from a highly back reef pool variable (~3-5°C of diurnal range) shows higher thermal tolerance than conspeci cs from an adjacent, less variable pool 6 . Rock Island reefs have a larger diurnal range (0.59°C) than the outer reefs (0.32°C; Fig. 1B), which could also facilitate thermal tolerance through plasticity mechanisms. In addition to Rock Island lineages being more thermally tolerant, outer reef-associated DB corals showed a trend toward fewer stress bands during the 1998 mass bleaching event when they had grown up and lived in the more variable Rock Islands (Fig.  5D). This suggests the variable Rock Island environments could also be promoting thermal tolerance via phenotypic plasticity. Alternatively, Rock Island conditions may have selected for the most thermally tolerant of the incoming DB larvae, or other Rock Islands conditions could mitigate bleaching in other ways. Additional studies that conduct reciprocal transplants of larvae and new recruits and robustly measure environmental variables would be needed to test hypotheses regarding the relative contributions of plasticity and adaptation in promoting thermal tolerance within the Rock Island habitats.
Future studies can also help determine the genetic underpinnings of the thermal tolerance differences seen across different lineages. For instance, Rock Island habitats differ from outer reefs in other parameters such as ow and light levels 11,16,37 . Shading in particular, has been posited as a potential explanation for the thermal tolerance of Rock Island coral communities as it may mitigate the irradiance stress that is often coupled with high temperature 11,38 . A lineage selected for lower light tolerance could appear thermally tolerant during heatwaves, as it may have experienced less stress residing in the shade. The LB lineage is most promising as a candidate to investigate a genetic basis for thermal tolerance in its own right because it maintains low bleaching levels across all habitats including unshaded outer reefs (Fig. 5). Adaptation to lower light levels could be a possible mechanism for the apparent thermal tolerance seen in the RD lineage, which is most common on one of the shadier Rock Island sites (Mecherchar; Fig. 2). Corals in Mecherchar grow mostly under the canopy of the trees along a narrow ledge (H. Rivera, personal observation). Of particular interest is the change in community composition seen between Mecherchar, which is RD dominated, and the nearby Mecherchar Channel site, which is LB dominated and where corals grow in the center of the pass, largely away from the shade of the nearby vegetation (Fig. 2). Understanding mechanisms of coral thermal tolerance would bene t from the kind of detailed, long term environmental data beyond temperature that is commonly available to terrestrial researchers but often lacking for marine environments.
Higher resolution of genetic and environmental data in coral systems will also help elucidate drivers of strong genetic divergence between sympatric lineages. Several recent studies suggest such patterns are the norm among coral species. For instance, along the Florida reef tract, both Siderastrea siderea and Montastraea cavernosa, show similar population structure to what we observe: highly diverged lineages, sometimes occurring sympatrically, and which differ across habitat types (in their case depth) 39 . In Panama, Orbicella faveolata populations harbor various distinct lineages, and these differ in thermal tolerance 40 . In Florida, Porites astreoides contains sympatric lineages that differ in thermal tolerance under experimental stress 7 . Acropora hyacinthus forms several strongly differentiated genetic lineages across mainland Japan and the Ryukyus archipelago, with one lineage appearing more adapted to colder temperatures and dominating the species' poleward range expansion in that region 41 . This recent work, along with our ndings, suggests that reef-building corals specialize to occupy narrow environmental niches, generating strong genetic differentiation even sympatrically and across small spatial scales. Whether these patterns may represent speciation in progress or simply be a characteristic of coral genetic diversity remains to be resolved and has important implications in the context of future species conservation efforts.
A central question within coral biology and conservation efforts is whether there are trade-offs between thermal tolerance and other key traits like growth or fecundity 32 . In addition to warm temperatures, the Rock Islands have pH and aragonite saturation levels near those expected in the open ocean in 2100 14 . Low pH and aragonite saturation can compromise coral growth and especially impacts skeletal density and facilitates bioerosion 27,28 . While other environmental conditions in the Rock Islands, such as light levels or turbidity can also in uence growth 42 , it is still worth examining potential growth and thermal tolerance trade-offs in the Rock Islands given their extreme pH and aragonite saturation conditions. LB corals were able to maintain high thermal tolerance and consistent growth regardless of pH conditions, indicating this lineage does not show any trade-off in its ability to handle multiple stressors (Fig. 5). In contrast, the outer reef-associated DB corals grew less, had lower density, and lower calci cation rates when found in the Rock Islands, where they show a trend toward higher thermal tolerance (Fig. 5). Though one could interpret the DB's lower growth as a trade-off with thermal tolerance, the challenging conditions for calci cation in the Rock Islands are more likely to be driving factors, especially since this lineage shows low thermal tolerance overall. The RD lineage, which is nearly exclusively found in the two lowest pH sites, Risong and Mercherchar 14,15 , shows lower growth metrics than the LB lineage, suggesting Rock Island conditions do have the potential to hinder coral growth (Figs. 1;3). Without being able to compare growth of the RD lineage under more favorable pH conditions, however, it is not possible to evaluate any trade-offs between its high thermal tolerance and growth. It appears that any combinations of the environmental factors across Rock Island and outer reefs sites do not affect the LB lineage in a substantial way, as it is able to maintain both growth and thermal tolerance across all habitats (Fig. 5). Thus, whatever trade-offs may exist, they do not appear to be ubiquitous.
The thermal tolerance of LB corals and their consistent growth across habitats have important implications for the future conservation and management of Palau's reefs and the use of thermally tolerant corals for reef restoration. Many reef systems are characterized by variable or warmer thermal regimes across space and through time, regimes that can select for and harbor thermally tolerant genotypes. Our results demonstrate that these environments can serve as breeding grounds for more tolerant corals (e.g. the LB and RD lineages) and that some of these (e.g. the LB lineage) can a) thrive and maintain their tolerance even when they disperse to cooler environments and b) maintain thermal tolerance without growth trade-offs. In addition, hydrodynamic models estimate that water exchange between Rock Island reefs and outer reefs began to slow only around 500 years ago 34 . The warmer temperatures and lower pH of the Rock Islands is in part due to long water residence times caused by limited ow 37 . Long residence times would also limit larval dispersal and increase selective pressures on local populations 34 . As such, it is possible that Rock Island corals adapted to warmer conditions in ~30-50 generations, assuming a generational time of 10-15 years for Porites cf. lobata 43 . This would imply that natural selection can increase coral thermal tolerance substantially over a much shorter time scale than normally thought, which could facilitate evolutionary rescue if such populations can disperse to more vulnerable areas.
As oceans worldwide continue to warm, corals derived from extreme habitats will be at a competitive advantage and may enable the survival of otherwise vulnerable reefs. Identifying and safeguarding natural breeding grounds of environmentally tolerant corals that can thrive under future climate conditions will be fundamental to the persistence of coral reef ecosystems worldwide in the coming decades. Nevertheless, the reality remains that curtailing climate-change and the greenhouse gas emissions that cause it will be the only way to truly safeguard our planet's biodiversity.
The site "Ngermid" has been referred to as "Nikko Bay" in previous publications. We use "Ngermid" here as that is the name preferred by Palauan natives. The "Outer Taoch" site contained samples from three outer reef locations: Airai (GPS coordinates: 7.33210, 134.56020, N= 5), Rael Dil (7.24990, 134.45073, N= 3), and a fringing reef (7.27193, 134.38115, N=30) immediately outside of Taoch Bay. The coordinates for this last site were used for mapping because most of the samples in this group are from this location. Due to the presence of various lineages within our dataset, population genetic metrics (e.g., F ST ) were not calculated by collection site, so this choice should have no bearing on results.

DNA extraction and Polymerase Chain Reaction (PCR) of microsatellites
Samples were processed as per 21 . A thawed ~1 mm 2 piece of coral was homogenized into a ne powder using a new standard safety razor blade sterilized with ethyl alcohol and amed. The homogenate was processed using the Qiagen® DNeasy Blood and Tissue DNA extraction kit according to the manufacturer instructions, with a modi ed Proteinase K incubation of at least 24 hours. Negative controls (N=5) without any coral tissue added were included every 70 samples and subjected to all the same downstream processing and analyses.

Microsatellite multi-locus genotyping
We used GENEMAPPER™ v3.0 (Applied Biosystems) to visualize electropherograms and call alleles. Scoring was conducted blind to the site of origin for each sample. The rst author scored alleles three separate times from scratch for all samples. Downstream analyses and results were consistent for all three sets. All automated allele calls were veri ed and curated manually to ensure accuracy and consistency between samples. After initial manual veri cation of all samples, raw allele sizes and allele call designations were exported and explored graphically. Boxplots of allele sizes by allele call designations were plotted for all markers, and all data points outside of the interquartile range were reveri ed manually and removed if peaks were of poor quality (i.e. very low height, non-standard shape, or possibly a spectral pull-up artifact from another channel that was not automatically detected). We also plotted allele size density curves by allele call designations, to identify samples with similar allele sizes but called as separate alleles. These techniques were used over several iterations to ensure allele calls were clean and of high quality.
For samples showing more than two alleles at a locus, the following steps were taken to select two alleles for population genetic analysis: I. For samples that were run more than once for any marker and showed a third allele in only one run, the singleton allele was dropped.
II. For samples run only once for a marker, the third allele was dropped if its height was less than half the second highest peak of the other two alleles.
III. For samples that were run more than once and showed 3 or more alleles consistently, or the sample was run only once but all alleles had roughly equal peak heights, the two alleles with the higher frequencies in the whole dataset were retained (i.e. the rarest allele(s) were dropped). This choice was made because it would be less likely to bias downstream analyses towards isolated populations.
Two microsatellite markers were dropped due to high rates of missing data (>40%). Samples with fewer than 10 of the remaining 12 loci were excluded (N=26); 322 samples were retained.

RAD-sequencing library preparation and sequencing
A subset of samples with su cient quality extracted DNA were then processed for RAD sequencing.
Genomic DNA concentrations were standardized using a Qubit™ 2.0 uorometer (Invitrogen) to 20 ng/ml. A total of 50 ml per sample was sent to Floragenex (Portland, Oregon) for single enzyme RAD library preparation with PstI enzyme digestion. Each sample was identi ed by a unique 10 nucleotide barcode. Samples (N=185) were sequenced as 100 base pair single end reads across 6 lanes of an Illumina Hiseq 4000™ using v4 chemistry at the University of Oregon Genomics Core facility.

RAD-seq processing and SNP-calling
Raw reads were processed using the 'process_radtags' module of Stacks v.1.46 46 , allowing for up to three mismatches in the sample barcode (this was the maximum number of mismatches at which the barcodes remained unique). Reads with low quality scores (PHRED<10) across a sliding window of 15% of the read length were discarded. We retained 78% of the original reads. Average sequencing depth was 8.5 million reads per sample. Three samples replicated within the plate showed 1.5-2-fold variability in sequencing depth. One sample which had an unusually high number of reads (>35 million) was discarded.
For read mapping and SNP calling, we used the dDocent pipeline 47 with a Porites lutea draft genome obtained from the REFUGE 2020 database (http://refuge2020.com/) as reference. dDocent clustered reads based on >95% similarity using CD-HIT 48 , mapped reads to the reference using the MEM algorithm of BWA 49 with a match score of 1, mismatch score of 3, and gap-open penalty of 4, and called SNPs using FreeBayes 50 with default values (E=3, m=PHRED 10, q=PHRED10, -V, and using the sampling sites as the populations designations). The resulting 'TotalRawSNPs.vcf' le was ltered using vcftools 51 and vc lter (https://github.com/jameshicks/vc lter) following the suggestions in the dDocent manual, with a nal thinning (-thin option in vcftools to keep only SNPs more than 150 bp apart, e.g. only one SNP per rad tag) to obtain a nal set of 12,761 bi-allelic SNPs in 146 retained samples.

Population structure
The R package adegenet 52 was used to explore genetic structure in both microsatellite and RAD-seq data. Principal component analyses were conducted using the dudi.pca() function on scaled and centered genind/genlight objects. We used the nd.clusters() functions to select the optimal number of groups based on Bayesian Information Criterion (BIC). Discriminant analyses of principal component (DAPC) was used using to visualize clusters, with the number of principal components retained determined through cross-validation xdapval() to avoid over tting. Nei's F ST was calculated using the gl.fst.pop() function from the package dartR 78 for RAD-seq data, using 100 bootstraps for estimating signi cance.
For microsatellite data, we used the function genet.dist() function from the package hierfstat 53 .
STRUCTURE v.2.3.4 54 was run using an admixture model with correlated allele frequencies and default parameters following previously used settings for corals 19,21 . Including sampling (geographic) information in the prior did not affect results and is not reported. MCMC chain settings were: 1 x 10 5 burnin, 1 x 10 6 iterations from K=1 to K=12, with 10 replicate chains per K.
We used the CLUMPAK feature 55 on the STRUCTURE Selector webserver to combine and visualize STRUCTURE output through the 'main pipeline' option with CLUMPP parameters: LargeKGreedy search, 10,000 random input orders, dynamic MCL, and default minimal cluster size. The Structure Selector webserver was used to run 'Best K' metrics which included methods to evaluate the optimal K described in 54,56,57 .

Temperature records and analyses
The Statistical metrics of temperature time series for each site were calculated using the 'zoo' 58 and 'xts' 59 packages in R. To examine daily temperature patterns, each time series was ltered in MATLAB 2015a using a bandpass Butterworth lter to retain signals between 5 and 30 hours in frequency and remove seasonal uctuations. Daily range (maximum-minimum temperatures) were calculated for each site in R.
Coral core sampling and analysis Coral cores were taken using an underwater pneumatic drill equipped with a diamond-tipped drill bit powered by compressed air from a SCUBA tank. Cores ranged from 10 to 204 cm long. Cores were dried in an oven and imaged using a Volume Zoom Helical Computerized Tomography (CT) Scanner at Woods Hole Oceanographic Institution. Scans were analyzed using an automated computer program developed and described in 28 and modi ed by 23 .
The presence/absence of stress bands during the 1998 and 2010 bleaching events (N=44, 5 sites), are data previously described and published in 13,28 , an additional 14 cores were analyzed for this study. Stress bands were de ned as a region of the core at least 1 mm thick in which density exceeded two standard deviations above the whole core average density, following the de nition in 13 .
While differences in stress band prevalence between RI and OR habitats had been previously shown 13 , we wished to test if genetic population groups provided additional explanation of these responses. A colony's genetic group was assigned as its predominant (>50%) STRUCTURE assigned group for K=4, which was the best K across several methods. Differences among lineages in growth parameters was tested using ANOVA and post-hoc Tukey tests. Differences in the proportion of cores showing stress bands between lineages was tested using a Chi-Squared test.