Selection of target sites

We aimed to collect coral specimens spanning coral phylogenetic diversity from a variety of Australian reefs. We targeted collection based on the 21 major coral clades defined in one of the most recent molecular phylogenies available at the start of the project33. Many of these monophyletic groups have since been defined as family-level taxa. Corals were collected at several sites on the east and west coast of Australia. These included Ningaloo Reef (Western Australia), Lizard Island, multiple reefs along the northern Great Barrier Reef, and Lorde Howe Island. Samples at Lizard Island were collected in both Summer and Winter, allowing for comparison of seasonal effects at one site across diverse corals.

Collection of metadata

During sampling, each coral, outgroup species, water, and sediment sample was associated with MIxS metadata42. This was accomplished by recording standardized metadata about each site prior to dives, and using an underwater metadata sheet (available at https://doi.org/10.6084/m9.figshare.5326870.v1). These metadata included basic features of coral species (as identified in the field), location, depth, water temperature, but also diver annotation of contact with macroalgae, turf algae or cyanobacteria (and the percent of the coral in contact); the presence of any visible tissue loss or disease signs; and coral color (using the Coral Reef Watch color charts43. Additionally, photographs of each coral were taken and released via openly accessible third party websites. They are easy to browse and thoroughly keyworded with taxonomy, location, and sample ID metadata on Flickr: https://flic.kr/s/aHsk9mjb54, and permanently archived in raw camera format with a spreadsheet linking filenames to colony names on FigShare at https://doi.org/10.6084/m9.figshare.5318236.v2.

Coral sampling

All coral samples were collected by AAUS-certified scientific divers, in accordance with local regulations. Relevant permit numbers are: CITES (PWS2014-AU-002155, 12US784243/9), Great Barrier Reef Marine Park Authority (G12/35236.1, G14/36788.1), Lord Howe Island Marine Park (LHIMP/R/2015/005), New South Wales Department of Primary Industries (P15/0072–1.0, OUT 15/11450), US Fish and Wildlife Service (2015LA1632527, 2015LA1703560), and Western Australia Department of Parks and Wildlife (SF010348, CE004874, ES002315). Only healthy corals were collected.

One goal of the project was to compare microbial diversity associated with the coral mucus, tissues and skeletons across many coral colonies. Each of these compartments represents a simplification of more complex structure, and much work remains to be done on the finer-scale distribution and dynamics of microorganisms across coral anatomy. For this project, we felt that a consistent reporting of these compartments across diverse corals represented a tractable step forward, given the scale of the project. Mucus was collected by gently agitating the surface of corals for ~30 s with a blunt 10 mL syringe. Exuded mucus or surface water (if no visible mucus was exuded) was then collected by suction. On the surface, settled mucus typically formed a distinct visible layer within the syringe. This was expelled into a cryogenic vial and stored in a dry shipper charged with liquid nitrogen for subsequent processing.

Tissue and skeletal samples were collected from each colony by hammer and chisel, or (for branching corals) by bone shears. These fragments were placed in sterile WhirlPaks and returned to the surface where they were snap frozen in a liquid nitrogen dry shipper until processing. In the laboratory, tissue was washed with sterile seawater (which removed visible mucus and detritus), then separated from skeleton using pressurized air of between approximately 800 and 2000 PSI (an ‘air gun’). Skeleton was sampled using a sterile chisel to isolate a ~1 cm3 region of skeleton that was not in direct contact with coral tissue. Skeleton samples were collected without regard to endolithic algae presence or absence (i.e., endolithic algae were neither specifically targeted nor excluded). Tissue slurries and skeleton samples were added directly to a MoBio PowerSoil Kit (MoBio Laboratories, Carlsbad, California) bead tube (which contains, among other things, a solution of guanidinium preservative) and stored at −80 °C until DNA extraction.

Sampling of reference samples

Because reef water and adjacent sediment might have an effect on the microbiota of corals from the same reef (especially in coral mucus), reef water and sediment were sampled at multiple sites. Surface seawater samples (1 L) were filtered through 0.22 μm Millipore Sterivex filters (Sigma-Aldrich, St. Louis, MO, USA) and reef sediment samples (2 mL) were collected in sterile cryogenic vials. Samples were snap frozen in a liquid nitrogen dry shipper, and subsequently stored at −80 °C until DNA extraction.

For comparison with corals from the same reef, we also opportunistically sampled non-scleractinian cnidarians from the genera Millepora (fire corals), Palythoa (zoanthids), Heliopora (blue corals), and Lobophytum (soft corals).

16S library preparation, sequencing, and initial quality control

DNA was extracted from skeleton, tissue, mucus, and environmental samples using the MoBio Powersoil DNA Isolation Kit. Two-stage amplicon PCR was performed on the V4 region of the 16S rRNA gene using the 515 F/806 R primer pair that targets bacterial and archaeal communities23. Extraction blank controls were also included in amplification and sequencing for quality assurance. The average concentration of extracted DNA used for PCR was 10.8 with a standard error of 1.2. First, 30 PCR cycles were performed using 515 F and 806 R primers (underlined) with linker sequences at the 5′ ends: 515F_link (5′-ACA CTG ACG ACA TGG TTC TAC A GT GCC AGC MGC CGC GGT AA −3′) and 806R_link (5′-TAC GGT AGC AGA GAC TTG GTC T GG ACT ACH VGG GTW TCT AAT −3′). Each 20 µL PCR reaction was prepared with 9 µL 5Prime HotMaster Mix (VWR International), 1 µL forward primer (10 µM), 1 µL reverse primer (10 µM), 1 µL template DNA, and 8 µL PCR-grade water. PCR amplifications consisted of a 3 min denaturation at 94 °C; 30 cycles of 45 s at 94 °C, 60 s at 50 °C and 90 s at 72 °C; and 10 min at 72 °C. Next, amplicons were barcoded with Fluidigm barcoded Illumina primers (8 cycles) and pooled in equal concentrations for sequencing. The amplicon pool was purified with AMPure XP beads and sequenced on the Illumina MiSeq sequencing platform (using V3 chemistry) at the DNA Services Facility at the University of Illinois at Chicago.

QIIME (v1.9)44 was used to process all 16S sequence libraries. Primer sequences were trimmed, paired-end reads merged, and QIIME’s default quality-control parameters used when splitting libraries. Chimeras were removed and 97%-similarity OTUs picked using USEARCH 7.045, QIIME’s subsampled open-reference OTU-picking protocol46, and the 97% Greengenes 13_8 reference database26. Taxonomy was assigned using UCLUST, and reads were aligned against the Greengenes database using PyNAST47. FastTreeMP48 was used to create a bacterial phylogeny with constraints defined by the Greengenes reference phylogeny. Following quality control, 9,441,738 usable reads remained. The number of per sample reads ranged from 2 to 38,523 with a median of 14,010, mean of 13,644, and standard deviation of 7565. Reads were partitioned across 129,305 unique OTUs (97% similarity cutoff). Sequencing success did not show any obvious trends with regards to host taxonomy or geographic location.

A ‘canonical’ rarefied OTU table was created and used for all downstream analyses except the linear model analyses. To create this table, OTUs were filtered out of the starting table if their representative sequences failed to align with PyNAST to the Greengenes database or if they were annotated as mitochondrial or chloroplast sequences. The beta_diversity_through_plots.py script was then used to rarefy the resulting table to exactly 1000 sequences per sample, and to calculate from this rarefied table multivariate dissimilarity measures including Bray-Curtis, Binary Jaccard, Weighted UniFrac, and Unweighted UniFrac. Also from this table, α-diversity statistics were calculated using alpha_rarefaction.py, including the number of OTUs observed, evenness, and Faith’s Phylogenetic Diversity.

Mitochondrial annotation and quality control

The primers used in this study were designed to selectively amplify the V4 region of bacterial and archaeal 16S rRNA gene, but we have noticed in many of our studies that they (and other standard primer sets) tend to strongly amplify corals’ mitochondrial 12S rRNA gene, which is the homolog of the bacterial 16S rRNA gene. Because our samples included species that were not used in the Huang and Roy 201520 phylogeny (including, critically, all outgroup taxa), these ‘off-target’ host mitochondrial reads were used to inform phylogenetic analyses and for an additional layer of quality-control. First, split_libraries_fastq.py was run on the raw forward reads without any quality trimming. Then, primers and adaptor sequences were removed, and USEARCH used to de-replicate 100% identical sequences. A frequency table was created and the data were filtered to contain only sequence variants with a total count of at least 100. Greengenes taxonomy was assigned to the remaining sequence variants with UCLUST as before, and sequence variants that had no match in the Greengenes database (e.g., putative non-bacterial or archaeal sequences) were isolated. For each host species, sequence variants were manually submitted to NCBI’s BLASTn web interface in order of their total abundance, comparing against the entire nr database. If a variant’s top 20 hits were annotated as coral mitochondria of any species, the sequence was copied to a FASTA file of host sequences. If all three compartments of a single coral individual of the same putative species still had unclassified variants that were more abundant, then manual annotation of those variants continued until either another coral mitochondrial sequence variant was found or there were no more variants in those samples that were more abundant than the previously annotated mitochondria. Using this method, no host sequences were found for some species of coral. The process was repeated for these species individually, without first discarding sequences that had counts of less than 100. In this way, mitochondrial sequences were eventually identified for every sample in the study.

Once all host species mitochondrial sequences were identified, the original frequency table of all unique sequence variants was filtered to contain only the identified host sequences. For each individual sample, the most abundant mitochondrial type was determined, and this information was then added to the sample’s metadata as its ‘12S genotype’. Then, all selected host sequences were aligned using MAFFT49 and de novo phylogenies were constructed in BEAST 2.4.250, with a chain length of 10 million, thinning interval of 1000, a log-normal relaxed clock model, and the site model selected using bModelTest51. The maximum clade credibility tree was selected using TreeAnnotator with a burn-in of 25% and common ancestor heights. This tree was compared to the expected topology (monophyletic Anthozoa, Hexacorallia, and Scleractinia, and otherwise matching the Huang and Roy 201520 molecular tree) to identify potential mismatches among the observed sequences and the field species identification. Regions of the tree with topology that differed from expectation were manually inspected.

Using this strategy, two coral individuals were noted whose field identifications placed them in the family Merulinidae, but whose sequence variants were strongly indicative of a relationship with the family Lobophylliidae. In these instances, further analyses verified that the same mitotype was detected in all three compartments of the same individual. Photos from collections in the field were consulted, and both were ultimately determined to have been misidentified in the field and in fact belonged to the genus Echinophyllia. Their metadata and annotations were updated to reflect this.

Aside from these two taxa, it was determined that unexpected topology in the de novo phylogeny was a result of imprecise resolution of the 12S V4 marker. For example, sequences from Millepora, Palythoa, and both octocoral species were placed in a monophyletic clade including the Complex corals, though they properly belong as outgroups to all Scleractinia. These errors emphasized the limitations of our opportunistic host sequence data to build a de novo phylogeny. Thus, having confirmed identifications of host species to within the resolution of the 12S marker, a new phylogeny was constructed with the topology constrained to exactly match the Huang and Roy 201520 molecular phylogeny and the known relationships of outgroup taxa. In cases where a single 12S genotype belonged to members of a polyphyletic group of taxa, we created separate tips for each monophyletic group. The mitochondrial sequence alignment and BEAST 2.4.2 were used to estimate relative branch lengths on this tree by supplying the starting tree and turning off all topology operators. The resulting tree was used for all phylogenetic analyses. As the branch lengths in this tree are derived from a relaxed clock model and limited sequence data, they are likely to represent some average between divergence times and degree of molecular evolution. Thus, analysis using these branch lengths represents a compromise between assuming correlation of traits is proportional to time since divergence and assuming that correlation of traits is proportional to overall evolutionary change since divergence.

Annotation of coral life history strategy

To assess connections between coral traits and microbiome structure, coral species sampled in this study were mapped to functional traits. These host features were added to the microbial mapping file, and used for tests of microbiome structure vs. host traits.

Coral life history strategies from Darling et al52. (‘weedy’, ‘competitive’, ‘stress-tolerant’, and ‘generalist’) were digitized and associated with coral species. Some species have recently been moved between genera based on updated phylogenetic evidence53. In these cases, both the original species name and the revised name are noted in the metadata. In some cases, species sampled were not annotated in Darling et al52. These were not assigned an annotation if annotated members of the same genus had mixed life-histories, or if only a single species of the same genus had been annotated. In cases where at least two members of the genus had been annotated and all annotated members shared the same life-history strategy, the same annotation was assigned to other members sampled from the genus.

Annotation of coral functional traits

Metadata associated with each species sampled was annotated with 28 reproductive, biogeographic, and morphological traits from the Coral Trait Database (CTDB) v. 1.1.123. These traits included basic details on coral distribution (abundance worldwide and on the Great Barrier Reef, range size, northerly and southerly limits, upper and lower depth limits), reproduction (sexual system, mode of larval development, propagule size, presence of Symbiodinium in propagules), phylogeny (genus and species ages), morphology (growth form, skeletal density, corallite maximum width, maximum growth rate) and conservation (IUCN Red List Category).

Adonis analysis of factors affecting microbial composition

We tested the influence of multiple host and environmental factors on the microbial community of each compartment. These results are presented in Fig. 1c and Supplementary Fig. 2, while the raw underlying data is presented in Supplementary Data S3. Throughout the analysis care was taken to account for the effects of rarefaction depth (we tested the robustness of the results at rarefaction depths of 1000 or 10,000 sequences/sample), β-diversity distance measure (we tested three distance measures), the degrees of freedom in each parameter (we used adjusted R2 values to account for differences in degrees of freedom), and to stringently control for the number of comparisons performed (using Bonferroni correction).

β-diversity distance matrices were calculated from separate OTU tables for coral mucus, tissue, and skeleton (outgroups and environmental samples were not included in this analysis). We calculated distance matrices using Weighted UniFrac distances, Unweighted UniFrac distances or Bray-Curtis dissimilarities. Then, for each host or environmental factor, the distance matrix was filtered to just those samples for which metadata were available (i.e., excluding ‘Unknown’ values). This prevented ‘Unknown’ values from being treated as a bona fide category in downstream statistical tests. The filtered distance matrix was then tested for clustering by factor using permutational tests (as implemented in Adonis in QIIME 1.9.1; 999 permutations per test). Because categories that can take on more values (e.g., species) are biased upwards in raw R2 values, we calculated adjusted R2 values for each category. These adjusted R2 values are primarily useful in that they allow for fair comparison between factors with differing degrees of freedom. Therefore, we present adjusted R2 values when comparing factors, but raw R2 values when discussing the percentage of variance explained (adjusted R2 values can no longer be interpreted as percent of variance explained). Importantly, we took care to separately filter the QIIME mapping file to exclude ‘Unknown’ values for each parameter under consideration. Failure to do so can lead to continuous variables (columns containing only numbers) being treated as categorical in QIIME, due to the presence of text values. This in turn can strongly influence inferred R2 values.

Summary of Adonis analysis of microbial beta-diversity

To present a summarized view of the Adonis analysis of microbial community beta-diversity, we compiled the R2 and p value obtained from each individual Adonis analysis. In Supplementary Fig. 1 the compiled adjusted R2 are presented in a heatmap. In Fig. 1c, we compared the relative influence of each host or environmental parameter on different host compartments by Z-score normalizing Adonis R2 values within columns. This has the effect of showing which compartments are most strongly influenced by a particular factor, independent of how influential that factor was overall. We present both views into the data because the unnormalized Adonis R2 values better emphasize the absolute magnitude of the microbiome response to each factor, whereas the Z-score normalized values better illustrate common patterns across compartments in host vs. environmental parameters. We emphasize that in both Fig. 1c. and Supplementary Fig. 2, clustering of rows and columns was performed without any prior specification of which factors were host vs. environmental. Thus, observed clustering of host vs. environmental factors emerges from features of the microbial communities themselves.

Machine learning analyses

All machine learning analyses were conducted through the supervised_classification.py script in QIIME (v1.9)44. This script implements random forest classification, which is a machine learning method for supervised classification. We used default parameters, which classify samples using inferred forests of 500 decision trees. We applied random forest classifiers to two tasks: 1) testing whether we can predict if a DNA sample came from coral mucus, tissue, or skeleton using microbial 16S rRNA data alone and 2) predicting whether within each coral compartment we can predict certain categorical features of a sampled coral using its microbiome (contact with turf algae, reef name, complex vs. robust clade membership, etc.). The results from random forest analysis of coral compartments are presented in the main text. The results for random forest analysis of host and environmental parameters are presented in Supplementary Data S6. Because error rates typically scale with both the number of categories—it is easier to predict the correct category for a binary category than one with 100 possibilities for example—we took care to consider the proportional increase in random classification relative to a baseline formed by random guessing. For both the compartment classification task and the trait classification task, we tested random forest classification on microbial phyla, orders or genera. For the compartment classification task we also tested random forest classification with 97% OTUs directly or predicted functional repertoires of coral microbiomes as inferred using the PICRUSt software. However, this was computationally expensive and yielded<0.1% improvement in classification error rates over classification based on microbial genera, and was therefore not pursued further.

Statistical analyses on the effect of phylogeny on the microbiome

Phylogenetic analyses were conducted in R v3.3.154. The packages ape (v3.5)55 and paleotree (v2.7)56 were used to manipulate trees and to calculate cophenetic distances. Univariate phylogenetic correlograms of α-diversity and distance-to-centroid measures were implemented using the package phylosignal (v1.1)57. Mantel tests and mantel correlograms of multivariate dissimilarities vs. phylogenetic distances were implemented using vegan (v2.5–2 ) 58. Size classes for the mantel correlograms were defined manually. Following Sturge’s rule, we created 11 distance classes. Due to the structure of the host phylogeny and the sampled species, four discrete phylogenetic distances greater than ~0.3 existed, corresponding to (1) all comparisons between the two major coral clades, (2) all comparisons between scleractinians and Palythoa (Zoantharia), (3) all comparisons between hexacorals and octocorals, and (4) all comparisons between anthozoans and Millepora (Hydrozoa). We first created distance classes that corresponded to each of these four discrete comparisons, then created the remaining seven distance classes by spacing them evenly across the smaller phylogenetic distances.

Phylogenetic Generalized Linear Mixed Models (pGLMMs) for analysis of the entire community were implemented using the package MCMC.OTU (1.0.10)59. MCMC.OTU wraps the package MCMCglmm (v2.24)60 to fit a model whereby sequencing depth is accounted for by using a sample’s total read count as the base level of a fixed per-OTU effect. Compositionality is accounted for by the inclusion of a per-sample random effect, and a number of other parameters and priors are set with defaults that are sensible for microbiome studies, such as ‘global effects’ of each specified factor (which control for and test effects on α-diversity) and independent error variance for each OTU. Analysis of the entire set of 97% OTUs was computationally impractical, so OTUs were first collapsed from the pre-rarefaction OTU table into their annotated genera using QIIME’s summarize_taxonomy.py. The package phyloseq (1.18.1)61 was then used to import and manipulate this table and its associated metadata. Samples with total counts less than 1000 or that were lacking relevant metadata were removed. The purgeOutliers command was applied to the data with an otu.cut value of 0.0001. For the first, more comprehensive GLMM, the command mcmc.otu was run with maximum corallite width, disease prevalence, and binary turf contact as fixed effects; geographic area, host phylogeny, and host identity as random effects; a chain length of 1,25,000, thinning interval of 5, and burn-in of 25,000; and with the inverse of the host phylogenetic covariance matrix supplied with the ginverse option. Subsequently, the command was run again with latitude and then coral colony size as the sole fixed effect, and with only host phylogeny as a random effect. Significance for each term was determined by calculating 95% credible intervals with HPDinterval and isolating those that did not include zero.

Cophylogenetic analyses

We reasoned that microbial groups that are most intimately associated with corals (whether commensal, mutualistic, or parasitic) are likely to have evolved in ways that led to patterns of cophylogeny with their hosts. A preliminary pipeline was developed to screen the microbiome for such groups. First, joined sequences were re-processed with the Minimum Entropy Decomposition (MED) pipeline34, discarding MED nodes with substantive abundances less than 100. Taxonomy was assigned to the resulting MED representative sequences as before with OTUs. Family-level groups of microbes were analyzed independently because higher taxonomic levels would be unlikely to have evolved within the same timescale as scleractinians, and lower taxonomic levels were more likely to contain misannotations. It was computationally impractical to analyze all microbial families, so only the most prevalent in each compartment were tested. An arbitrary threshold of 50% prevalence was chosen. For each family in each compartment, all MED nodes were isolated that were the most abundant representative in at least one sample. This conservative approach was done partly out of concern that spurious sequences generated by sequencing error could influence the downstream phylogenetic analyses, and partly to reduce each dataset to a size that was practical for phylogenetic inference and GLMMs. The representative sequences were then combined from these nodes with reference sequences for each family. Reference sequences were randomly subsampled from the Greengenes 13_8 99% OTU database such that each dataset contained the MED nodes of interest, 75 random full-length 16Ssequences belonging to the family of interest, plus 10 random ‘outgroup’ sequences belonging to any other family from the same order.

Each collection of sequences was then aligned using MAFFT in QIIME. Phylogenetic trees were built using BEAST 2.4.2 with a chain length of 100 million, thinning interval of 1000, a log-normal relaxed clock model, and the site model selected using bModelTest. The maximum clade credibility tree for each group was selected using TreeAnnotator with a burn-in of 25% and common ancestor heights.

A separate pGLMM was then fit for each microbial family in each tissue compartment. The raw MED table was imported into R using phyloseq and filtered to contain only samples with counts greater than 1000. The resulting table was merged with each microbial family’s phylogenetic tree using phyloseq, a process that automatically filters all sequences from the table that are not represented on the tree. Samples were further filtered from this table if they did not retain a count of least 10. Phylogenetic covariance matrices based on the bacterial and host phylogenies were then generated62. Phylogenetic covariance matrices based on the bacterial and host phylogenies were generated using the function inverseA on each host tree. The Kronecker product of the resulting matrices was then computed for use as the ‘coevolutionary’ covariance matrix. The Kronecker product of each phylogenetic covariance matrix and an identity matrix was computed for use as microbial identity x host phylogeny and microbial phylogeny × host identity interaction effects62

Binary models were fit with MCMCglmm using a single fixed effect of the log of the sequencing depth, ‘global’ random effects of host phylogeny, host identity, microbial phylogeny, and microbial identity, all combinations of host-by-microbe phylogenetic and identity random interaction effects, and a geographic area-by-microbial identity random interaction effect. Altogether this approach is similar to the models described in reference62. Our models were fit with a chain length of 1,250,000, thinning interval of 50, and burn-in of 250,000. After the model was fit, convergence was assessed by verifying that the Effective Sample Sizes (ESS) of all covariance terms were greater than 200. Intraclass correlation coefficients (ICCs) were calculated for each iteration, with 95% credible intervals calculated with HPDinterval. Factors with ICC lower credible bounds greater than 0.01 were considered significant.

To independently analyze subclades of Endozoicomonas-like bacteria, a custom QIIME-formatted taxonomy database was created with sequence annotations corresponding to clades C, HS, HS-R, and HS-C from the initial analysis. Taxonomy was then assigned to all MED nodes and Endozoicimonaceae Greengenes reference sequences using UCLUST with max-accepts set to 1. The above procedure of filtering, selecting reference sequences, building a phylogeny, and fitting pGLMMs, was then repeated based on each annotated subclade instead of each family. HS-R within only Robust clade corals was also analyzed by first filtering other samples from the dataset.

Code availability

Analysis code is available on GitHub: https://github.com/zaneveld/GCMP_Australia_Coevolution