Diverse coral reef invertebrates exhibit patterns of phylosymbiosis.

Microbiome assemblages of plants and animals often show a degree of correlation with host phylogeny; an eco-evolutionary pattern known as phylosymbiosis. Using 16S rRNA gene sequencing to profile the microbiome, paired with COI, 18S rRNA and ITS1 host phylogenies, phylosymbiosis was investigated in four groups of coral reef invertebrates (scleractinian corals, octocorals, sponges and ascidians). We tested three commonly used metrics to evaluate the extent of phylosymbiosis: (a) intraspecific versus interspecific microbiome variation, (b) topological comparisons between host phylogeny and hierarchical clustering (dendrogram) of host-associated microbial communities, and (c) correlation of host phylogenetic distance with microbial community dissimilarity. In all instances, intraspecific variation in microbiome composition was significantly lower than interspecific variation. Similarly, topological congruency between host phylogeny and the associated microbial dendrogram was more significant than would be expected by chance across all groups, except when using unweighted UniFrac distance (compared with weighted UniFrac and Bray-Curtis dissimilarity). Interestingly, all but the ascidians showed a significant positive correlation between host phylogenetic distance and associated microbial dissimilarity. Our findings provide new perspectives on the diverse nature of marine phylosymbioses and the complex roles of the microbiome in the evolution of marine invertebrates.

Introduction

Phylosymbiosis occurs when microbial community relationships reflect the evolutionary history of the host [1-3]. The term was first coined to describe the impact of a host phylogenetic signal on gut microbial community relationships in Nasonia parasitoid wasps [2, 4], and the phenomenon has since been investigated in a diverse range of taxa and environments, e.g., the gut microbiomes of mammals and insects [1, 5, 6], the skin microbiome of ungulates [7], the endolithic microbiome of coral [8] and the root microbiome of plants [9]. These studies have confirmed that phylosymbiosis occurs in the simplest as well as the most diverse microbial communities and the discovery of virus/host phylosymbioses [10] demonstrates that the phenomenon is not limited to prokaryotes. As phylosymbiosis has become more frequently observed, the mechanisms underpinning these patterns are of increasing interest.

Evolutionary processes such as codivergence and coevolution are distinct from phylosymbiosis, establishing the need of an alternative term [1]. Namely, phylosymbiosis is a pattern observed at one moment in time and space, which does not assume a stable evolutionary association between a host and its microbiota or congruent ancestral splits, nor does it assume vertical transmission of microbial symbionts [11]. While it is possible that different evolutionary processes contribute to the mechanisms behind phylosymbiosis [8, 12], complex and dynamic systems that acquire high numbers of microbes from the environment are likely structured by other mechanisms. For example, horizontal transmission of microbes filtered through phylogenetically congruent host traits, biogeography of a host and the microbiota, and dispersal of microbes among conspecifics all potentially contribute to observed phylosymbiosis patterns [12-15]. These explanations are not necessarily mutually exclusive. Within a complex microbiome where both vertical and horizontal transmission occurs among obligate and facultative microbial members, phylosymbiosis is expected to rely on multiple mechanisms [3, 6].

Despite the extensive literature supporting phylosymbiotic relationships, host phylogeny does not always correlate with microbial community (dis)similarity. For example, in contrast to other mammals, no significant congruence was observed between skin microbiome composition and host phylogeny in the case of carnivores [7]. Similarly, no phylosymbiotic signal could be detected in the case of the intestinal microbiota of 59 Neotropical birds [16] and the gut microbiomes of bats are more similar to birds than other mammals [17]. There are multiple reasons why phylosymbiosis may not occur. First, factors such as environment and diet may obscure phylosymbiotic signals, which have been successfully controlled for in some studies [1, 4]. Second, in some cases, host genotype exerts strong effects on microbiome composition that are independent of host phylogeny [8, 18, 19]. Finally, host physiology can structure the microbiome [20], however, physiological traits may not always be consistent with host phylogeny [21]. Therefore, patterns of phylosymbiosis may be dependent on a certain host taxonomic level (i.e., host family), where host genotype effects are reduced and host physiological traits and phylogeny are congruent.

Reef invertebrates provide interesting opportunities for testing hypotheses of phylosymbiosis, as they often host diverse microbial communities acquired by combinations of vertical and horizontal transmission [22-25] that can be dynamic among different environments [26, 27]. Here, we first characterise the microbiomes of four groups of coral reef invertebrates: scleractinian corals, octocorals, sponges and ascidians. We then test three recommended analyses to investigate phylosymbiosis: (a) comparison of intraspecific and interspecific variation in microbiome composition, (b) comparison of the topology of host phylogeny and hierarchical clustering of its associated microbial community, and (c) correlation of host phylogenetic distance with microbial community dissimilarity [3, 14]. We hypothesise that a phylosymbiotic signal will be found across all four groups to show that host phylogeny is a dominant factor in microbiome structure of reef invertebrates. Through an improved understanding of microbial community dynamics using phylosymbiosis, our knowledge of how a microbiome is structured and maintained in complex marine holobionts will be enhanced [25].

Materials and methods

Sample collection

Tissue samples from 3 to 5 replicates of 30 species of coral reef invertebrates (12 corals, 10 octocorals, 5 sponges and 3 ascidians) were collected on SCUBA from seven locations across the central and northern sectors of the Great Barrier Reef (GBR) (Table S1; Fig. S1). On sampling trips to Broadhurst Reef, Davies Reef and Orpheus Island, August 2017 (Table S1), adult colonies no larger than 30 × 30 cm were collected using hammer and chisel and returned to the reef after sampling. Alternatively, sampling of invertebrates was performed in situ. On the surface, colonies/samples were isolated and placed in running seawater (0–2 h) until processing. Each invertebrate was sampled for 3–5 fragments ~5 cm in length using either a hammer and chisel or dive knife (coral), or sterile razor blades (all other invertebrates). In addition, seawater samples were collected from the central GBR sites in August 2017 as an environmental control (Table S1). All samples were collected under the marine parks permits G12/35236.1 and G15/37574.1

Sample processing and preservation

Fragments were rinsed in autoclaved calcium- and magnesium-free seawater (CMFSW; NaCl: 26.2 g, KCl: 0.75 g, Na2SO4: 1 g, NaHCO3: 0.042 g, per 1 L) to remove any loosely attached microbes. For scleractinian coral, tissue was removed from the skeleton by pressurised air into ~30 ml of CMFSW. Coral blastate was homogenised by vortex for 1 min and 2 × 2 ml aliquots were kept for DNA extraction. Aliquots were centrifuged for 10 min at 10,000 × g, the supernatant was removed, and tissue pellet was either snap frozen in liquid nitrogen or preserved in 1 ml dimethyl sulfoxide-EDTA salt saturated solution (DESS) and kept at −80 °C (Table S1). For octocorals and sponges, fragments were cut into small pieces ~0.5 × 0.5 cm3 using a sterile razor blade, snap frozen in a 2 ml cryovial and stored at −80 °C until DNA extraction. Alternatively, a 15 ml falcon tube with ~7 ml DESS was filled with the dissected tissue until approximately a 1:1 ratio of tissue:DESS was reached. The ascidians Lissoclinum patella and Polycarpa aurata were dissected longitudinally and the tunic layer was removed and snap frozen as described above. Colonies of the remaining ascidian Didemnum molle were dissected into three equal parts as the tunic was too small to isolate and preserved in 1 ml DESS and kept at −80 °C. Seawater was collected from each site (excluding the Ribbon Reefs (RR) and Osprey Reef) ~1 m above the benthos at the area of sample collection using 4 × 5 L retractable water bottles (washed and sterilised with 10% hydrochloric acid). Approximately 2–3 L were then filtered through 0.22 µm sterivex filters and stored at −80 °C (where −80 °C was not available, samples were stored at −20 °C for 1–5 days before being transferred to −80 °C upon returning to the lab).

DNA extraction and sequencing

Approximately 0.05 g of tissue was used for DNA extraction using the DNeasy PowerBiofilm Kit (QIAGEN Pty Ltd, VIC Australia 3148). Extraction was performed following the manufacturers protocol with the BioSpec Mini-Beadbeater-96 used for mechanical lysis at 3–5 cycles of 30–60 s depending on the difficulty to break down the tissue. Genomic DNA was sent to the Ramaciotti Centre for Genomics (UNSW, Sydney Australia) for 16S rRNA amplicon sequencing on the Illumina MiSeq platform using the modified V4 region primer set, 515F (GTGYCAGCMGCCGCGGTAA) [28] and 806R (GGACTACNVGGGTWTCAAT) [29]. Samples were prepared for sequencing with the Earth Microbiome Project’s 16S Illumina Amplicon protocol and sequencing was performed following the standard Illumina protocol for 16S rRNA gene amplicon library prep. Sequencing of the host phylogenetic markers COI, 18S rRNA and ITS1 was performed at the Beijing Genome Institute following the BGISEQ-500 library prep protocol on the BGISEQ-SE400 module. COI (~712 bp), 18 S (~470 bp) and ITS1 (~288 bp) were amplified using the primer pairs, LCO1490 (GGTCAACAAATCATAAAGATATTGG) and HCO2198 (TAAACTTCAGGGTGACCAAAAAATCA) for COI [30] and V4_18S_Next.For (CCAGCASCYGCGGTAATTCC) and V4_18S_Next.Rev.B (ACTBTCGYTCTTGATYARNGA) were modified from Pirredda et al. [31] for 18S rRNA. For ITS1, the custom primers 18S-F1759 (GGTGAACCTGCGGAWGGATC) and 5.8S-R40 (CGCASYTDGCTGCGTTCTTC) were designed by retrieving all available sequences from our target species and aligning them using MAFFT [32]. Full length barcodes were assembled from single-end 400 bp reads using the HIFI-SE pipeline [33].

16S rRNA gene amplicon analysis

Sequences were analysed using QIIME2 (v 2018.4) [34] by first demultiplexing reads and denoising following the DADA2 pipeline [35]. Taxonomic assignment was performed using a Naive Bayes classifier pre-trained on the Silva 132 99% OTU database modified to the V4 region primer set 515F/806R. The resulting amplicon sequence variant (ASV) table was filtered for chloroplast, mitochondrial and eukaryotic sequences. A phylogenetic tree was reconstructed using the qiime fragment-insertion sepp command (QIIME2 v 2019.1), which places the ASVs into a larger, well-curated 16S rRNA reference phylogeny containing >200,000 representative tips (GreenGenes 13.8, 99% OTU) [36]. The resulting tree was then trimmed to the original reference sequences and used for subsequent UniFrac analyses. ASV and taxonomic tables were imported into R studio v.3.5.0 [37] for further analysis with extensive use of the packages ‘phyloseq’ [38], ‘vegan’ [39], ‘ggplot2’ [40], ‘ggtree’ [41], ‘ape’ [42], ‘phangorn’ [43] and ‘dplyr’ [44].

Characterisation of microbial diversity and composition

The following analyses were conducted at the ASV level, excluding visual representations of relative abundance. Relative abundance for each microbial phylum was calculated and grouped by invertebrate taxonomy to give a broad overview of microbial profiles of each invertebrate group. In addition, the top 25 most abundant microbial families across the entire dataset were shown to give an overview of the lower taxonomic levels. As the taxonomic profile of the blanks was sufficiently different from the marine invertebrate profiles, with only 0.4% of sequences present in the top 25 family level ASVs, these samples were removed from further analysis. Rarefaction curves were calculated and plotted to illustrate the total diversity of ASVs captured against the sampling effort. Alpha diversity was calculated using both species richness (total number of ASVs retrieved per sample) and Shannon–Wiener diversity index on a dataset rarefied to 3500 sequences (equal to the sample with the lowest number of sequences). Beta diversity was calculated on non-rarefied data using the Bray–Curtis dissimilarity measure by first standardising the data by the species maximum and then by the sample total (Wisconsin double standardisation). This method of normalisation was chosen for beta diversity as transforming data to proportions returns the most accurate Bray–Curtis dissimilarities [45]. The resulting dissimilarity scores were visualised using non-metric multidimensional scaling (NMDS) to observe overall patterns in microbial community structure among the different invertebrates. Analysis of variance (ANOVA) and a post-hoc Tukey’s test with unplanned comparisons and a Bonferroni correction were used for significance testing of alpha diversity, while permutational multivariate analysis of variance (PERMANOVA) was used for beta diversity using the pairwiseAdonis function for post-hoc analysis.

Host phylogenetic reconstructions

Representative sequences for COI, 18S rRNA and ITS1 from each species in each taxonomic group were aligned separately using MUSCLE [46] and then concatenated using DAMBE [47]. Concatenated octocoral and sponge alignments were further curated using Gblocks [48] to remove poorly aligned, high gap regions. Evolutionary model selection was performed using JModelTest2 (Supplementary Table 2) [49] and phylogenetic analysis was conducted in Mr Bayes v3.2.7 [50] using the outgroups Carteriospongia foliascens for corals, octocorals and ascidians and Cladiella sp. for sponges. Outgroups were selected based on their low phylogenetic relatedness to the ingroup and low variability in microbiome composition among sample replicates. Evolutionary history was inferred using Bayesian inference with the Markov Chain Monte Carlo method using two independent runs of 5,000,000 generations and all models converged at <0.01.

Phylosymbiosis analysis

The 16S rRNA gene dataset was subsampled to each taxonomic group and analysed independently. Gorgonians did not contain enough species within our dataset to compare host phylogeny with microbial composition and were added to the soft coral dataset to create an octocoral group. Intraspecific against interspecific variability of microbiome composition was compared using pairwise comparisons of Bray–Curtis dissimilarity between each sample. Welches t test was used for significance testing following an arcsine transformation to normalise the 0–1 distribution, while an ANOVA and post-hoc Tukey’s test with unplanned comparisons and a Bonferroni correction were used to test for significant differences in intraspecific variation among invertebrate groups.

Microbial dendrograms were built in QIIME2 using the qiime diversity beta-rarefaction command. Within each invertebrate ASV table subset, all ASVs that appear two times or less and those that are present in only one sample were removed to reduce noise from potentially spurious and transient ASVs. Each sample was then pooled by host species and rarefied over 1000 iterations to the host species with the lowest number of reads following the method of Brooks et al. [1]. Hierarchical clustering of host species from the resulting table was performed using the UPGMA clustering method based on Bray–Curtis dissimilarity and both weighted and unweighted UniFrac distances. Microbial dendrograms along with phylogenetic trees and pooled ASV tables were imported into R studio for analysis.

To assess topological congruency, host phylogenetic tree topology was compared with the microbial dendrograms using the normalised Robinson–Foulds (nRF) metric, where 0 is complete congruence and 1 is no congruence. Branch lengths were removed in host phylogenetic trees for visualisation and a significance value was calculated using the RFmeasures function [14] with 9999 permutations. Correlation between host phylogenetic distance and microbial dissimilarity was analysed by first creating a distance matrix of pairwise phylogenetic distances between each host species and distance matrices of Bray–Curtis dissimilarity and weighted and unweighted UniFrac distances using the pooled ASV tables. A Mantel test was used to test for correlation between host and microbial distance matrices using Pearson correlation with 9999 permutations. A similarity percentages (SIMPER) analysis was used to identify which ASVs were contributing to dissimilarity between host species that showed incongruence.

Results

Sample collection and sequencing

Field collections resulted in a total of 161 samples across 30 species of reef invertebrates (Table S1). In addition, eight seawater samples, two blank extractions and two sequencing positive controls were sequenced. For 16S rRNA amplicon sequencing, this yielded a total of 10,415,183 reads in 173 samples, which was reduced to 8,611,147 high quality reads following quality control and denoising. For host phylogeny, successful COI sequences were obtained for all 30 species, however, 18S rRNA sequencing was unsuccessful for Acropora formosa, Acropora hyacinthus, Diploastrea heliopora, Heteroxenia sp. and Isis hippuris and ITS1 sequencing was unsuccessful for Lissoclinum patella and Didemnum mole. As a result, ITS1 was not used for ascidian phylogeny.

Rarefaction curves for each sample approached asymptotes, illustrating that total ASV richness for each sample was captured (Fig. S2). However, rarefaction to the sample with the lowest number of reads (Isis hippuris: 3323 reads; excluding blanks) resulted in a loss in diversity in some samples. Nonetheless, overall trends showed that both ASV richness and ASV diversity (Shannon–Wiener Index) were both significantly different across the broad taxonomic associations (richness; ANOVA; F(5, 163) = 7.01, p < 0.001; Fig. 1) (Shannon diversity; ANOVA; F(5, 163) = 4.64, p < 0.001; Fig. 1). Post-hoc comparisons revealed that seawater had a significantly higher ASV richness than the ascidians (p = 0.024), while coral had a significantly higher ASV richness than ascidians (p = 0.006), soft corals (p = 0.006) and sponges (p = 0.003). For ASV diversity, post-hoc comparisons revealed an increase in diversity in coral compared with the ascidians (p = 0.009) and soft corals (p = 0.046), and an increase in seawater compared with the ascidians (p = 0.014). However, unlike richness, no difference was seen in ASV diversity between corals and sponges (p = 1.0).

Fig. 1

ASV richness (top panel) and Shannon–Wiener diversity index (bottom panel) for each invertebrate group and seawater.

Letters indicate groups which are significantly different from each other.

A total of 62 microbial phyla were observed across the invertebrate groups and microbial profiles showed a high degree of uniformity at the phylum level. Microbial taxonomy mentioned here and herein are ASV sequences affiliated to that taxonomic classification, with Proteobacteria, Cyanobacteria and Bacteroidetes among the dominant phyla across all marine invertebrates (Fig. S3). However, differences were evident even at the broad taxonomic level, with the octocorals (soft coral and gorgonians) hosting a higher relative abundance of Tenericutes (mean = 4.71% ± 1.63 SE and 11.12% ± 6.28 SE, respectively) compared with other invertebrates, while sponges were associated with more Chloroflexi, Acidobacteria and Cyanobacteria (mean = 19.09% ± 2.29 SE, 9.86% ± 1.64 SE and 28.31% ± 3.67 SE, respectively).

Relative abundance at the family level indicated far more variation in taxonomic profiles among the invertebrate groups (Fig. 2). The three groups of anthozoans (coral, soft coral and gorgonian) were clearly different to the other marine invertebrate classifications and mostly dominated by the common Endozoicomonadaceae (mean = 33.52% ± 4.19 SE, 38.41% ± 4.58 SE and 42.88% ± 10.74 SE, respectively). Sponges consisted of a high relative abundance of Cyanobiaceae (mean = 27.87% ± 3.74 SE), comprised of the commonly found cyanobacteria Prochlorococcus and Synechococcus (Silva database classification), as did seawater (mean = 32.26% ± 2.46 SE). Ascidians appeared more variable, with Rhodobacteraceae, Porticoccaceae, Cyclobacteriaceae and unclassified Alphaproteobacteria, all abundant within the top 25 bacteria at the family level.

Fig. 2

Relative abundance of the top 25 prokaryotic families found across each invertebrate group as well as seawater and blank extractions.

Bubble size is proportional to the relative abundance of each prokaryotic family (y-axis) within a host group (x-axis).

Between sample variability (beta diversity) showed there was an overall weak clustering of samples by their broad taxonomic classifications (Fig. 3). Particularly the three anthozoans (coral, soft coral and gorgonian) and ascidians had low homogeneity in microbial composition. Comparatively, sponge and seawater samples formed clusters that indicated consistent microbial composition across samples. Microbial composition was confirmed statistically to be associated with host taxonomy (PERMANOVA; F(5, 163) = 2.58, p < 0.001), however, only a small amount of variation in the data was explained by the broad taxonomic classification (R2 = 0.073). When samples were instead grouped by host species, the amount of variation explained increased dramatically (PERMANOVA; F(30, 138) = 2.01, R2 = 0.30, p < 0.001). Lastly, beta-diversity analysis showed there was a significant association to collection site (PERMANOVA; F(6, 162) = 1.90, R2 = 0.066, p < 0.001), however, only a small amount of variation could be explained by this variable, and since many species were collected from only one reef, it is likely the variation is due to species-specific microbiomes.

Fig. 3

Bray–Curtis dissimilarity based on microbial composition visualised using NMDS.

Each symbol represents a sample where colour is the associated host and shape is reef zone where sample was collected.

Assessment of phylosymbiosis among coral reef invertebrates

All four marine invertebrate groups showed lower intraspecific Bray–Curtis dissimilarity in microbial composition compared with interspecific Bray–Curtis dissimilarity (coral: t(364) = 13.53, p < 0.001; octocoral: t(302) = 18.84, p < 0.001; sponge: t(200) = 34.80, p < 0.001; ascidian: t(69) = 19.09, p < 0.001), confirming lower microbiome variability among conspecifics (Fig. 4). Furthermore, intraspecific variation was significantly different among the invertebrate groups (ANOVA; F(3, 818) = 231.15, p < 0.001), with the exception of the ascidians and octocorals (t = 1.85, p = 0.39), highlighting sponges and coral with the highest and lowest microbiome homogeneity, respectively.

Fig. 4

Intraspecific and interspecific Bray–Curtis dissimilarity scores for each invertebrate group.

Interspecific variation (red boxplots) in the microbiome was significantly greater than intraspecific variaion (blue boxplots) for each invertebrate group.

Comparing the topology of host phylogenetic trees with the corresponding microbial dendrograms (nRF test) and measuring the correlation of host phylogenetic distance with microbial dissimilarity (Mantel test) further revealed significant levels of phylosymbiosis across all four groups of invertebrates (Table 1). Patterns of phylosymbiosis were significant in sponges using all tests and metrics (Figs. 5a and S4), while Bray–Curtis and weighted UniFrac metrics found significant patterns of phylosymbiosis using the nRF and Mantel tests in corals (Figs. 5b and S5a) and octocorals (Figs. 5c and S6a). Using the unweighted UniFrac distance, phylosymbiosis patterns were significant only using the Mantel test but not the nRF test for coral (Fig. S5b) and octocoral (Fig. S6b) and no patterns were detected in the ascidians (Fig. S7b). Perfect congruency between host phylogeny and microbial dendrograms was observed in the ascidians using both the Bray–Curtis and weighted UniFrac metrics (Figs. 5d and S7a). Despite this, no significant phylosymbiosis was observed using the Mantel test. This opposing result is likely due to the low sample size combined with marked differences in microbial composition among the three ascidians (Fig. S8).

Table 1

Normalised Robinson–Foulds (nRF) and mantel statistics across Bray–Curtis, weighted and unweighted UniFrac beta-diversity metrics.

	Bray–Curtis	Weighted UF	Unweighted UF
Sponge
nRF	RF = 0.02, p < 0.001	RF = 0.4, p = 0.006	RF = 0.4, p = 0.01
Mantel	r = 0.71, p < 0.001	r = 0.78, p = 0.006	r = 0.75, p = 0.03
Coral
nRF	RF = 0.69, p < 0.001	RF = 0.69, p < 0.001	RF = 0.92, p = 0.15
Mantel	r = 0.37, p = 0.02	r = 0.38, p = 0.01	r = 0.42, p = 0.03
Octocoral
nRF	RF = 0.64, p < 0.001	rRF = 0.82, p = 0.02	RF = 0.91, p = 0.24
Mantel	r = 0.23, p < 0.001	r = 0.36, p < 0.001	r = 0.25, p < 0.001
Ascidian
nRF	RF = 0, p < 0.001	RF = 0, p < 0.001	RF = 0.5, p = 0.34
Mantel	r = −0.03, p = 0.63	r = 0.46, p = 0.17	r = 0.18, p = 0.46

Fig. 5

Host phylogeny and microbial dendrogram comparisons for each invertebrate group.

a–d Host phylogenies are inferred from COI, 18S rRNA and ITS1 sequences while microbial dendrograms are based on Bray–Curtis dissimilarity for microbial composition of each host species. Cladiella sp. was used as an outgroup for a sponges, while C. foliascens was used as an outgroup for b coral, c octocoral and d ascidians. Numbers at nodes reflect posterior probability for clade support in host trees and jackknife support values in dendrograms. Branch tips are coloured to reflect clades in host phylogeny. Initials in brackets next to species names refer to collection site. BR Broadhurst Reef, DR Davies Reef, OR Osprey Reef, PI Palm Islands (Orpheus and Pelorus), PR Pandora Reef, RB Rib Reef, RR Ribbon Reefs. P. massive refers to massive Porites sp.

A select few species were collected from multiple locations and showed contrasting results in relation to phylosymbiosis. The sponge Ircinia ramosa and octocoral Sarcophyton sp. were collected from two locations and both correctly formed a clade with their conspecifics (Figs. 5, S4 and S6), which was supported by uniform microbial profiles (Figs. S9 and S10). Conversely, the octocoral Sinularia sp. and the coral species Porites cylindrica and Seriatopora hystrix did not form clades with their conspecifics from different locations and there was a reduced overall phylosymbiotic signal (Figs. 5, S5 and S6). A SIMPER analysis revealed that shifts in the relative abundance of ASVs assigned to Endozoicomonadaceae were consistently the top contributors to the dissimilarity observed between species collected from two sites (Table S3; Fig. S11). For example, Porites cyclindrica collected from the Palm Islands (PI) had a dramatic reduction in Endozoicomonadaceae compared with those collected from the RR, where the mean relative abundance of Endozoicomonadaceae fell from 82.9% (±4.32 SE) to 3.31% (±1.69 SE). Similarly, the microbial profile of Sinularia collected from RR differed from the two Sinularia species collected from PI, with colonies from RR hosting a lower relative abundance of Endozoicomonadaceae and a higher relative abundance of unknown bacteria and Fusobacteriaceae (Fig. S11).

Additional incongruences were observed among the groups where sample location was not a factor. The overwhelming majority of extant corals fall into one of two major clades, the Robusta and Complexa. This split was only partially reflected in the Bray–Curtis and weighted UniFrac microbial dendrograms, although in most cases, species within a genus or family clustered together (Figs. 5b and S5). Similarly, host phylogeny was recapitulated in the microbiome of only certain clades of octocorals using Bray–Curtis and weighted UniFrac metrics, such as the microbiome of Briareum and species within the family Alcyoniidae (Sarcophyton, Sinularia and Cladiella), with the exception of Sinularia collected from the RR (Figs. 5c and S6). However, no congruence was seen between gorgonian phylogeny and microbial composition, which can again be attributed to ASVs assigned to Endozoicomonadaceae (Table S3; Fig. S10). Lastly, although the signal of phylosymbiosis in sponges was strong and robust across all analyses, the main incongruence was due to an unclassified Ircinia sp., which did not form a clade with its sister species in the host phylogeny (Fig. 5a), and highlights the unresolved phylogenetic relationships among the Ircinia [51].

Discussion

This study evaluates the signal of phylosymbiosis in diverse coral reef invertebrates, finding evidence that host evolutionary history helps shape the microbiome in sponges, corals, octocorals and ascidians. By testing three commonly used methods for phylosymbiosis analysis, we show that all groups have lower intraspecies microbiome variability compared with interspecies. This was combined with greater topological congruency between host phylogeny and the microbial dendrogram than would be expected by chance, except when using the unweighted UniFrac distance in corals, octocorals and ascidians. Interestingly, all invertebrate groups but the ascidians exhibited a significant correlation between host phylogenetic distance and microbial dissimilarity across all beta-diversity metrics.

Our results demonstrate that sponges have a strong signature of phylosymbiosis, which likely reflects the uniform microbiome structure in sponges compared with other coral reef invertebrates [18]. This was observed through low intraspecific variation and high homogeneity in the microbiome when the same species was collected from different reefs. Sponges are also known to have a relatively stable microbiome in response to temporal variation and environmental perturbations [52-54]. A stable microbiome may lead to a strong phylosymbiotic signal if there is less influence from the surrounding environment, leaving host factors to be the primary structuring element of the microbiome [53]. Importantly, while sister species were included in the analysis, overall the sponges sampled here span a larger phylogenetic diversity compared with the other groups, which may increase the chance to observe phylosymbiosis. Our results agree with previous conclusions of a significant correlation between host phylogeny and microbiome dissimilarity and validate a prominent role of host phylogeny in shaping the sponge microbiome [18, 55].

A signal of phylosymbiosis was demonstrated in coral, which was characterised by a tendency of corals of the same genus or family to cluster together. However, incongruences were observed where the same species was collected from two different locations, primarily due to a shift in the relative abundance of Endozoicomonadaceae. Shifts in Endozoicomonadaceae have been documented previously, normally in response to host stress [26, 56]. As shifts in the microbial community can often precede visual signs of an unhealthy holobiont [57, 58], it is plausible the decrease in Endozoicomonadaceae is linked to an unknown event. Second, coral tissue samples are often contaminated by the coral mucus, which is known to have a dynamic microbial community shifting in composition between new and aged mucus [59]. However, bacteria within the tissues of corals are housed within coral-associated microbial aggregates and these communities likely have a more stable association with the host [60, 61]. Therefore, developing approaches to target tissue-specific microbes could be beneficial to understanding phylosymbiosis and other questions related to microbial symbiosis in corals.

Similar clustering of coral microbiomes has been observed in Caribbean corals. This partially reflected coral phylogeny, as congenerics showed comparatively low microbial dissimilarity and the two major coral clades tended to cluster together, however, inconsistencies were seen when looking at the species level [62], and reflect the results seen here on the GBR. Further evidence of phylosymbiosis in coral was found in an analysis of 691 coral samples collected Australia wide [8]. The endolithic microbial community showed the strongest signal and was the best predictor of the deep phylogeny between the Robusta and Complexa clades. Tissue microbiomes also illustrated evidence of phylosymbiosis, however, the signal was absent in the coral’s surface mucus layer. This emphasises an increasing strength of phylosymbiosis where direct environmental factors are reduced. In addition, a small number of microbial lineages, including those within Endozoicomonadaceae, demonstrated co-phylogeny with their host, while other clades had a more generalist host distribution. It is possible that host-specialist clades play a minor role in phylosymbiosis through codivergence and future work should aim to untangle the mechanisms behind phylosymbiosis [14].

Research on the microbiome structure of octocorals is limited compared with corals, and we show for the first-time direct evidence for phylosymbiosis. The phylosymbiotic signal in octocorals was similar to corals and incongruences also occurred when there was a shift in the relative abundance of Endozoicomonadaceae. Octocorals are known to have a more stable and less diverse microbial community than hard corals [63], consistent with our finding that overall microbial diversity was lower and microbiome uniformity higher in octocorals compared with hard corals. While this likely influences the phylosymbiotic signal, a direct comparison between octocorals and corals (and other invertebrate groups) cannot be drawn due to the differences in phylogenetic relatedness between host species. Furthermore, the phylogenetic markers used in this study were chosen to capture both mitochondrial and nuclear evolution across a broad range of diverse species. However, octocorals have poorly understood phylogenetic relationships, with little concordance between morphological, nuclear and mitochondrial data [64]. The incorporation of alternative phylogenetic markers optimised for each taxonomic group may further improve the analyses of phylosymbiosis and comparisons among groups. Finally, octocoral identification in the field is extremely challenging especially when trying to resolve to species level [65]. Despite these limitations, we still observe a significant signal of phylosymbiosis, which is likely to strengthen with improved phylogenetic relationships and species identification.

Ascidians showed complete congruence between the host phylogeny and microbial dendrogram for both Bray–Curtis and weighted UniFrac metrics, yet no correlation existed between host phylogenetic distance and microbial dissimilarity. Our results therefore do not provide strong support for phylosymbiosis in the group, yet they highlight the need for multiple lines of evidence when evaluating phylosymbiosis [3]. For example, we find that when sample numbers are low, particularly when marked changes are observed among the microbiomes of host species, the dendrogram approach was more sensitive to patterns of phylosymbiosis compared with the Mantel test. Furthermore, unweighted UniFrac methods were unable to identify a phylosymbiotic signal in the ascidians and had the least power to identify a signal across all invertebrate groups, which agrees with previous conclusions on weighted and unweighted beta-diversity metrics [14]. As this method does not account for the abundance of ASVs, it is less likely to identify beta-diversity patterns in highly diverse microbiomes that are dominated by a relatively small number of bacteria

Our study overwhelmingly found that host phylogeny is reflected in the microbiome of marine invertebrates, particularly notable when considering several confounding factors. Sampling of the reef invertebrates occurred over four field trips that spanned a 1-year timeframe, potentially obscuring phylosymbiosis patterns due to seasonal influences on the microbiome [66]. Furthermore, our samples are from wild colonies collected from multiple locations on the GBR which introduces local environmental differences including water quality and the pelagic communities that serve as host diet. Preservation methods also varied across organisms, including snap freezing and the use of salt saturated dimethyl sulfoxide-EDTA. While these preservation approaches have been shown to have little effect on the microbial composition of coral, it could have influenced alpha diversity [67]. Finally, sample representation differed among the four groups and likely has an important impact on the strength of the phylosymbiosis signal. For example, only three ascidian species (and one outgroup) were used, whereas four related species are recommended [3]. Had more species been included in the analysis, with a larger number of taxonomic sister species, a more reliable representation of phylosymbiosis would likely have been achieved.

This is the first study to systematically assess phylosymbiosis among diverse groups of marine invertebrates. We identified a phylosymbiotic signal across all invertebrate groups with multiple methods, of which sponges consistently showed a significant signal using all beta-diversity metrics. Increased intraspecific variability of the microbiome in both scleractinian corals and octocorals was often associated with a change in the relative abundance of Endozoicomonadaceae. This microbial family is characterised by host-specialist and host-generalist clades and is assumed to be a dynamic member of the coral holobiont [8]. Host-specialist clades may contribute to phylosymbiosis in corals and octocorals through codivergence, while host-generalist clades obscure the signal through host infidelity. Here, we provide a foundation to begin exploring the mechanisms behind phylosymbiosis and further our understanding on host-microbe symbiosis and coevolution in marine invertebrates.

Supplementary tables

Supplementary figures