Coral sampling
Between 2011-18, we collected tissue from 543 Porites cf. lobata colonies using a hammer and chisel while on SCUBA. Tissue was preserved in RNAlater™ (Invitrogen, Waltham, MA), incubated overnight at 4°C, and frozen at -20°C (N=329), or frozen directly at -80°C (N=20), or preserved in 95% ethanol and frozen at -20°C (N=194) until DNA extraction. Colonies were sampled haphazardly within each site, across 13 sites (Fig. 2), based on morphological characteristics of Porites lobata detailed in 44.
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., FST) 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 per21. A thawed ~1 mm2 piece of coral was homogenized into a fine powder using a new standard safety razor blade sterilized with ethyl alcohol and flamed. The homogenate was processed using the Qiagen® DNeasy Blood and Tissue DNA extraction kit according to the manufacturer instructions, with a modified 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.
We amplified 14 microsatellite markers with fluorescently labeled primers developed for P. cf. lobata by21,45. PCR settings were: (1) initial denaturation at 94°C for 5 min, (2) 35 cycles of [94°C for 20 seconds, annealing at 52, 54, or 56°C (plex-dependent) for 20 seconds, 72°C for 30 seconds], and (3) final extension for 30 minutes at 72°C. The Pennsylvania State University Nucleic Acid Facility measured fragments on an ABI 3730 (GeneScan) with a LIZ-500 internal size standard.
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 first 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 verified and curated manually to ensure accuracy and consistency between samples. After initial manual verification 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 re-verified 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:
- 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.
- 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.
- 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 sufficient quality extracted DNA were then processed for RAD sequencing. Genomic DNA concentrations were standardized using a Qubit™ 2.0 fluorometer (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 identified 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.4646, 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 pipeline47 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-HIT48, mapped reads to the reference using the MEM algorithm of BWA49 with a match score of 1, mismatch score of 3, and gap-open penalty of 4, and called SNPs using FreeBayes50 with default values (E=3, m=PHRED 10, q=PHRED10, -V, and using the sampling sites as the populations designations). The resulting ‘TotalRawSNPs.vcf’ file was filtered using vcftools 51 and vcffilter (https://github.com/jameshicks/vcffilter) following the suggestions in the dDocent manual, with a final 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 final set of 12,761 bi-allelic SNPs in 146 retained samples.
Population structure
The R package adegenet52 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 find.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 overfitting. Nei’s FST was calculated using the gl.fst.pop() function from the package dartR78 for RAD-seq data, using 100 bootstraps for estimating significance. For microsatellite data, we used the function genet.dist() function from the package hierfstat53.
STRUCTURE v.2.3.454 was run using an admixture model with correlated allele frequencies and default parameters following previously used settings for corals19,21. Including sampling (geographic) information in the prior did not affect results and is not reported. MCMC chain settings were: 1 x 105 burn-in, 1 x 106 iterations from K=1 to K=12, with 10 replicate chains per K.
We used the CLUMPAK feature55 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 in54,56,57.
Temperature records and analyses
The Coral Reef Research Foundation (http://wtc.coralreefpalau.org/) provided 30-minute interval in situ temperature data for Mecherchar, Helen, Drop Off, Ngerdiluches, Ngerchelong, Ngermid, and Kayangel. These data were recorded using U22 loggers (Onset Technologies, MA). Temperature data covered periods from 2010-2017 and from 2-15 meters depth. The foundation indicated accuracy was determined to be within 0.1°C through pre- and post-deployment calibration against a NIST traceable mercury thermometer and that individual thermographs were also cross calibrated with each other. Temperatures for Risong and Taoch were obtained from U22 loggers (Onset Technologies, MA) deployed between 2-5 meters depth by the Cohen Lab at Woods Hole Oceanographic Institution from 2011-2013 recording at 15-minute intervals.
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 filtered in MATLAB 2015a using a bandpass Butterworth filter to retain signals between 5 and 30 hours in frequency and remove seasonal fluctuations. 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 in28 and modified by23.
The presence/absence of stress bands during the 1998 and 2010 bleaching events (N=44, 5 sites), are data previously described and published in13,28, an additional 14 cores were analyzed for this study. Stress bands were defined 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 definition in13.
While differences in stress band prevalence between RI and OR habitats had been previously shown13, 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.
Symbiodiniaceae inter-transcribed spacer-2 (ITS2) denaturing gradient gel electrophoresis (DGGE) genotyping and sequencing
We analyzed 27 coral samples representing all four lineages, across 8 outer reef and Rock Island sites. To test whether the dominant symbiont community of coral colonies shifted across Palau’s strong environmental gradients, we amplified the ITS2 region of Symbiodiniaceae’s nuclear ribosomal DNA and visualized bands using DGGE following protocols in60,61. Briefly, the ‘ITSintfor2’ and ‘ITS2clamp’ primers were used for initial amplification with a touchdown PCR protocol consisting of: (1) initial denaturation at 94°C for 2 minutes, (2) 20 cycles of [94°C for 20 seconds, initial annealing temperature of 62°C for 10 seconds and decreasing at 0.5°C intervals every cycle until 52°C, then 68°C for 30 seconds], (3) continuing with another 18 cycles at annealing temperature of 52°C, and (4) a final extension for 10 minutes at 68°C.
Products were loaded onto 8% acrylamide gels with a 40-75% denaturing gradient and run for 24 hours at 90 volts. The gels were stained in 1 liter of deionized water with 10 ml of SYBR Red™ (ThermoFisher) for 30 minutes, de-stained in 1 liter of deionized water for 30 minutes, and then visualized with a UV gel imager. DNA from Cladocopium (C15) cultures was obtained from the LaJeunesse Laboratory (Pennsylvania State University) and run alongside Palauan samples for band identification.
Representatives of any additional bands seen in Palauan samples were excised using a sterile pipette tip, homogenized in 5ml of molecular grade water, and reamplified using the ‘ITS2intfor2’ and ‘ITS2rev’ primers and a standard PCR protocol: (1) initial denaturation at 92°C for 3 minutes, (2) 35 cycles of [92°C for 30 seconds, annealing at 52°C for 40 seconds, 72°C for 30 seconds], and (3) final extension for 10 minutes at 72°C. Products were visualized on a 1% agarose TAE gel, and successfully re-amplified samples were purified using a MinElute™ PCR Cleanup Kit (Qiagen) and sent for Sanger sequencing at Sequegen (Worchester, MA). Sequences were then aligned to Symbiodiniaceae sequences on the NCBI ‘nt’ database using the MEGABLAST algorithm with default parameters on the NCBI BLAST webserver.