Algal isolates of three species of Symbiodinium, formerly known as clade A, were used for this study (Table 1). Cultures were grown in ASP-8A media [pH 8.5, 35ppt, Blank (1987)] at control (26°C) and high (32°C) temperature conditions. Cultures were grown on Percival incubators and placed on a glass shelf on top of full-spectrum fluorescent lights at 100 μmol quanta m^–2 s^–1 facing up, with a 12:12 light/dark cycle photoperiod. Light intensities were measured inside a 300 ml culture flask with media using a 4π sensor (Biospherical, United States). Cultures were maintained by adjusting cell density to 1 × 10⁵ cells ml^–1 every 10 days by adding fresh media, acclimated to the control temperature (26°C) for at least 2 months (Supplementary Figure 1A). To measure cell density, a sterile rubber policeman was used to scrape all cells from the flask walls to ensure collection of all cells. The culture was homogenized manually and 1 ml aliquot was collected into a 2 ml screwcap tube. To detach cell aggregations, lugol was added (200 μl) and tubes were vortexed for 5 min. Cells were then counted manually with a hemocytometer. To begin the experiment, replicate 15 ml clear glass tubes containing 8 ml of culture with 1 × 10⁵ cells ml^–1 were placed vertically on top of a glass shelf with full spectrum lights facing up, in an incubator set at 26 or 32°C (n = 5 per isolate per temperature) (Supplementary Figure 1B). Tubes were undisturbed until after 27 days of incubation, when F_v/F_m was measured for all replicates (n = 5 per isolate per temperature) and cells were collected for DNA extraction (n = 3 per isolate per temperature). Cells were collected by scraping all cells attached to the glass tube with a disposable inoculating loop and centrifuged at approximately 2,000 g for 1 min, removing the supernatant and flash freezing with liquid nitrogen immediately. DNA was isolated using the MO BIO PowerSoil DNA isolation Kit following manufacturer’s instructions. A previous study established that S. microadiaticum (strain CassKB8) and S. necroappetens (strain RT80) have intermediate thermotolerance, since they were able to grow at high temperature but not as well as at normal temperature. S. pilosum (strain RT130) grew at the same rate in normal and high temperature conditions, suggesting that it has higher thermotolerance (Díaz-Almeyda et al., 2017). S. pilosum and S. necroappetens were collected in the Caribbean and are part of the same isolate collection (formerly known as Trench collection) while S. microadriaticum was isolated in the North Pacific (formerly known as BURR collection) (Table 1).

Samples were amplified using specific primers pairs (Caporaso et al., 2011) for the V4 region of the 16S rRNA gene. A 25 μL PCR reaction was prepared with 10 μL 5Prime HotMaster Mix (QuantaBio), 0.5 μL forward primer 515F (GTGCCAGCMGCCGCGCGGTAA) (10 μM), 0.5 μL reverse primer 806R (GGACTACHVGGGTWTCTAAT) (10 μM), 1 μL of BSA (10 mg/ml), 1 μL template DNA (10 ng), and 12 μ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. Samples were submitted for sequencing to the Department of Energy’s Joint Genome Institute (JGI) and were processed according to their iTag sequencing protocol. Individual libraries were prepared and barcoded before being multiplexed and sequenced on an Illumina MiSeq sequencer machine in a 2 × 300 run mode.

Sequences were demultiplexed and processed according to the JGI’s iTag amplicon sequencing protocol for taxonomic identification. Initial processing of 16S rRNA gene sequence libraries was performed using the Quantitative Insights Into Microbial Ecology (QIIME2; 2019.7.0) package (Bolyen et al., 2019). Read trimming, filtering, denoising, and amplicon sequence variant (ASV) inference was performed with the DADA2 implementation of QIIME2 (Callahan et al., 2016). Taxonomy was assigned to ASVs using the Naive Bayes classifier pre-trained on the Silva database version 138 clustered at 99% sequence similarity and optimized for the 515F/806R primer pairs. Extraction blanks and template PCR controls were not included as part of our sequencing effort. We therefore applied additional filters to the reads to remove possible contaminants. Initially, ASVs were removed if they were not present in a minimum of three of the 18 samples sequenced. Additionally, low abundance ASVs totaling less than 63 reads (0.2 percent of the total reads from the sample with the fewest reads) were removed from each sample. The additional filtering removed rare ASVs, many of which were present in a single sample and highly likely to be contamination of the culture. Filtering reduced the total ASV count from 255 to 124. Phylogenetic diversity metrics were generated with alignments performed with MAFFT and trees generated with FastTree (Supplementary File 1), which was followed by alpha and beta diversity calculations.

Differential abundance testing was performed with DESeq2 [v1.32.0, Love et al. (2014)], against sample abundances normalized by the geometric mean of positive counts. Fold change values were adjusted for dispersion using the adaptive shrinkage estimator from the ashr package [v2.2-47, Stephens et al. (2022)], and ASVs differentially abundant beyond a false discovery rate of 0.01 were collected (Supplementary Table 1). Contributions of each ASV to Bray-Curtis dissimilarity were calculated for using the SIMPER function from the vegan package (v2.5-7) (Supplementary Table 2). Permutational multivariate analysis of variance (PERMANOVA, 999 permutations) was performed as well (Supplementary Table 3). Initially, samples were rarefied to 30,000 reads and pairwise comparisons were performed between 26 and 32°C (Supplementary Table 4).

Data visualization, including PCoA plots of unweighted and weighted UniFrac distance results calculated in QIIME2 were generated with qiime2R (v0.99.6)¹ and phyloseq [v1.36.0, McMurdie and Holmes (2013)] R packages (RStudio Team, 2020). Overlapping taxa (grouped at the genus level) between species and temperature were identified using the online tool from VIB/Ugent Bioinformatics and Evolutionary Genomics Website.²

The microbial communities of each Symbiodinium sp. were characterized via sequencing of the V4 region of the 16S rRNA gene. We recovered an average of 50,623 raw reads per sample. After filtering and processing the sequences, the average number of sequences per sample was 48,507, with a minimum of 31,090 and maximum of 68,716. After removal of low abundance taxa, we identified 124 bacterial ASVs across all samples and we summarized them at the genus level by replicate (Supplementary Figure 2) and by treatment and species (Supplementary Figure 3). Samples grouped largely by Symbiodinium species when clustered by Bray-Curtis dissimilarity (Figure 1). Permutational multivariate analysis of variance (PERMANOVA) confirmed diversity composition significantly varied by species (999 permutations, q-value < 0.05), indicating Symbiodinium species harbor specific microbial communities (Supplementary Table 3).

Comparison of community composition with SIMPER identified less than 20 ASVs driving 85% of the observed dissimilarity between species (Supplementary Table 2). Additionally, OM182_clade, which accounted for the greatest percentage of dissimilarity between S. microadriaticum and S. necroappetens, was absent from S. microadriaticum, while present in both S. necroappetens and S. pilosum (Figure 2). Labrenzia accounted for 18.9% of dissimilarity between S. pilosum and both other Symbiodinium species. These findings were further corroborated by differential abundance analysis with DESeq2. We found the aforementioned taxa to be differentially abundant between Symbiodinium species (Supplementary Table 1).

A closer inspection of these ASVs revealed that they belong to a small subset of the genera found across the three species, including Muricauda, Balneola, Labrenzia, Phaeodactylibacter, Sphingorhabdus, and OM182_clade, each with >8% influence in respective species comparisons (Figure 3 and Supplementary Table 2). Notably, four genera (Muricauda, Labrenzia, Roseovarius, and the uncultured SM1A02 group of the family Phycisphaeraceae) were commonly found across the three Symbiodinium species at both 26 and 32°C, with each dinoflagellate taxa associating with presumably a unique strain represented by a single ASV.

Symbiodinium microadriaticum harbored the greatest number of unique microbial genera (37 at 26°C and 36 at 32°C), whereas S. pilosum harbored the least, with seven species specific genera at both 26 and 32°C (Figure 4). No additional genera were shared between S. pilosum and S. necroappetens beyond the aforementioned four taxa, whereas S. microadriaticum and S. necroappetens shared 13 to 14 genera at 26 and 32°C, respectively.

While the four shared genera were often found to exhibit high relative abundance within the samples, other taxa were represented in greater abundance within each Symbiodinium species. In S. microadriaticum, Phaeodactylibacter (15.21–16.75%) and usually a member of the OM60(NOR5) clade of Halieaceae (6.08–15.80%) were the dominant genera. In S. necroappetens, OM182_clade of Gammaproteobacteria (17.51–21.96%), Balneola (13.89–24.78%), Marinobacter (5.89–17.18%), and Labrenzia (11.14–18.01%) were the dominant genera. In S. pilosum, Muricauda (15.09–21.61%), Sphingorhabdus (8.60–18.06%), and Labrenzia (38.32–49.74%) were the dominant genera (Figure 3).

Shannon Diversity Index (H’) values were calculated for each sample (Figure 5). S. microadriaticum had the highest average H’ with no significant differences between temperatures (4.67 at 26°C and 4.78 at 32°C), followed by S. necroappetens (H’ = 3.60 at 26°C and 3.61 at 32°C). S. pilosum had the lowest average H’ of 3.04 at 26°C and 2.67 at 32°C. One-way ANOVA analysis showed significant differences between species [F_(5,12) = 218.3432, p < 0.0001] but only S. pilosum had a significant change in H’ between the two temperatures (Tukey HSD, p = 0.0053).

Our samples did not cluster well by temperature (q-value > 0.05). We did, however, find ASVs that were significantly differentially abundant between 26 and 32°C. Seventeen ASVs were differentially abundant in cultures of S. microadriaticum. Two of the differentially abundant taxa belonged to the shared genera Labrenzia and Muricauda, and increased by 8- and 2-fold (log₂), respectively. In S. necroappetens, both bacterial genera exhibited differential abundance, with Labrenzia increasing 5-fold (log₂) and Muricauda decreasing in relative abundance by 5-fold (log₂). Nine ASVs were differentially abundant in S. pilosum, although none belonged to the four genera commonly found between the three dinoflagellate species.

All three Symbiodinium species shared members of Muricauda (Flavobacteriaceae), Labrenzia (Rhodobacteraceae), Roseovarius (Rhodobacteraceae), and a member of Phycisphaeraceae (SM1A02). These findings are consistent with previous studies of Symbiodiniaceae-associated microbial communities (Lawson et al., 2018; Camp et al., 2020). Importantly, although microbial communities grouped by Symbiodinium species, they did not group phylogenetically by Symbiodinium species. Rather, they grouped by location and collection of origin (Supplementary Figure 4). S. necroappetens and S. microadriaticum are more closely related phylogenetically than S. pilosum (LaJeunesse et al., 2015). However, S. necroappetens had a microbial community more similar to S. pilosum than to S. microadriaticum. S. necroappetens was isolated by Robert K. Trench at least 20 years ago (LaJeunesse, 2001), S. pilosum was isolated at least 27 years ago (Iglesias-Prieto and Trench, 1994), and both were isolated from the Caribbean by the same laboratory with protocols which included adding antibiotics and using a low nutrient culture media to remove bacteria (Schoenberg and Trench, 1980). This may explain the low bacterial diversity we observed in S. pilosum. In contrast, S. microadriaticum was isolated at least 35 years ago by a different laboratory using f/2 media (Taguchi and Kinzie, 2001). It is therefore unknown whether grouping of bacterial communities is due to long-term culturing artifacts or true symbiotic associations, as these differing protocols may have resulted in very different founder effects, leading to the differences we observed. Maire et al. (2021) found in situ samples of Symbiodiniaceae differed greatly from corresponding cultures, suggesting that co-culturing may have influenced our results rather than true symbiosis. The same study also found the presence of intracellular bacteria in Symbiodiniaceae to be common, diverse, and highly distinct from extracellular bacteria. We are unable to separately characterize the responses of intracellular and extracellular communities to thermal stress because we did not distinguish these communities. With these considerations, our results provide insight into the potential connection between Symbiodinium thermotolerance and associated bacterial communities, and point toward future studies on in situ Symbiodinium that distinguish between closely- and loosely-associated bacteria.

Bacterial diversity of S. microadriaticum and S. necroappetens remained stable at high temperatures while that of S. pilosum was significantly lower at 32°C. Stable H’ values across different environmental conditions can indicate a lack of substantive change in community composition, or that the combined richness and evenness of the community is relatively stable (Kimbro and Grosholz, 2006; Willis, 2019). The flexibility that we observed in S. pilosum may be related to its higher thermotolerance, as acclimating to rapid changes in environmental conditions has been associated with dramatic shifts in bacterial community composition (Vargas et al., 2020; Voolstra and Ziegler, 2020). In contrast to our results, Camp et al. (2020) found that the microbiome of Durusdinium trenchii, a thermotolerant Symbiodiniaceae, remained stable at high temperatures.

The bacterial communities of S. microadriaticum and S. pilosum also clustered by temperature. S. necroappetens samples clustered together but not separately by temperature, most likely due to experiencing a much smaller change with temperature. These results support the idea that microbiome responses to heat stress are unique to Symbiodinium species, perhaps regardless of dinoflagellate thermotolerance. Our results also suggest that increased thermotolerance may correlate with lower levels of microbial diversity and greater individual variability in community composition. According to the “Anna Karenina principle,” the increased prevalence of environmental stressors may reduce a host’s ability to regulate its microbiome composition (Zaneveld et al., 2017). While our results did not offer evidence for this principle, if well-supported with future research, the Anna Karenina principle has strong implications for coral reefs in the context of climate change, as a strong link between rising ocean temperatures and abundance of coral pathogens within communities of coral-associated microbes has been described (Tout et al., 2015), and microbial communities without functional redundancy can be subject to environmental stressors such as ocean warming (Voolstra and Ziegler, 2020). Further research using controlled studies to investigate this relationship may find a link between the thermotolerance of a microbial community and instability of the holobiont under high temperatures. The potential benefits of greater flexibility and diversity in microbiome community composition for thermotolerance should be further explored in future studies using larger sample sizes and other Symbiodinium species. Additionally, while thermotolerance is known to be affected by other factors such as genetics, environmental factors, and host-symbiont type mismatches (Brown et al., 2002; Mieog et al., 2009), the degree of contribution (positive or negative) by certain microbial taxa, their potential symbiotic interactions with Symbiodinium, and their relationships to Symbiodinium thermotolerance remain unexplored and should be considered for future studies.

No ASVs consistently increased with temperature across all Symbiodinium species. Rather, different ASVs were differentially abundant by temperature for each Symbiodinium species. While some or all of these differences may be a result of long-term culturing, we offer insights on the potential importance of key ASVs for Symbiodinium thermotolerance.

We found Labrenzia to be particularly abundant across all Symbiodinium species, but especially in S. pilosum, the most thermotolerant species. Additionally, Labrenzia were not differentially abundant between temperature treatments in S. pilosum but were significantly less abundant at 32°C in the two intermediately thermotolerant Symbiodinium. Labrenzia are known for their ability to produce dimethylsulfoniopropionate, which scavenges reactive oxygen species and thus may play a role in increasing Symbiodinium stress tolerance (Diaz et al., 2016; Curson et al., 2017). It is also plausible that these bacteria could be directly reducing pathogenic microbes, which might opportunistically take advantage of heat stressed Symbiodinium. An oyster-colonizing Labrenzia sp. exhibited broad antimicrobial activity and reduced the growth of a pathogenic Roseovarius (Amiri Moghaddam et al., 2018), and labrenzbactin, a siderophore produced by a species of Labrenzia, has been shown to exhibit antimicrobial activity in vitro (Sharma et al., 2019). However, its effect on Symbiodinium-specific pathogens has not yet been explored. The increased prevalence of Labrenzia in S. pilosum and its relatively stable abundance with temperature stress may be related to thermotolerance, or it may be a result of having lower levels of diversity in general. Research on the ecological and physiological roles of Labrenzia for Symbiodinium fitness is necessary to discern whether these bacteria are generally mutualistic, parasitic, or commensal partners in heat-stressful conditions.

Roseovarius were also present in all three host species, in much lower relative abundances than Labrenzia, and was highest in S. necroappetens. Other studies have similarly found Roseovarius in relatively low abundance with Symbioiniaceae (Camp et al., 2020; Maire et al., 2021). Notably, Pootakham et al. (2019) found that Roseovarius were one of the two most prevalent genera associated with heat-stressed corals. Members of Roseovarius have also been implicated in coral diseases such as black band disease (Ng et al., 2015). However, we did not find that Roseovarius significantly increased in relative abundance with increased temperature in Symbiodinium (Figure 2). Our results therefore suggest that Roseovarius may not play a similar role in heat stressed Symbiodinium, or perhaps only in more thermosensitive species. Further research on the functional ecology of Roseovarius with Symbiodinium under heat stress is necessary to investigate these hypotheses.

Muricauda was also shared by all three Symbiodinium species we cultured and were differentially abundant by temperature in the two intermediately thermotolerant Symbiodinium, but not in S. pilosum. Several studies have demonstrated that Muricauda are intimately associated with Symbiodiniaceae (Han et al., 2016; Lawson et al., 2018). In the first comparison of intracellular, closely-associated, and loosely-associated bacterial communities, Maire et al. (2021) found Muricauda present in most in situ samples of different Symbiodiniaceae species, but consistently in less than 1% abundance. We found Muricauda in greater than 1% abundance in all Symbiodinium, but especially in S. pilosum, in which Muricauda consistently accounted for over 15% of each sample. This discrepancy may be the result of culturing, as cultured samples have previously been found to differ substantially from in situ Symbiodiniaceae samples (Maire et al., 2021). Regardless of its potentially low abundance in wild Symbiodiniaceae, Muricauda is one of 12 bacterial genera ubiquitously associated with various reef-building corals across the globe, suggesting its functional importance not only to Symbiodiniaceae but also to the coral holobiont (Bernasconi et al., 2019). Furthermore, Muricauda may serve a functional role in Symbiodinium heat resistance. A close relative to Muricauda lutaonensis enhanced thermal and light stress resistance of cultured Symbiodiniaceae (Motone et al., 2020). We did not find significant differences in Muricauda by temperature for any Symbiodinium species, but in situ sequencing may indicate whether they play a key role in Symbiodinium thermotolerance.

Sphinghorabdus were one of the dominant taxa in S. pilosum but were not present in greater than 1% abundance in the intermediately thermotolerant host species (Figure 3). Sphinghorabdus have been isolated from Gorgonian coral and have been successfully cultivated (Keller-Costa et al., 2017; Silva et al., 2018), but otherwise are not well-known, so further research may reveal their importance for Symbiodinium thermotolerance.

In S. pilosum, a member of Phycisphaeraceae (SM1A02) accounted for the greatest amount of community dissimilarity between temperature treatments, displaying greater abundance at 26°C. This may be due to strong competitive pressure from the highly dominant Labrenzia in thermally-stressed S. pilosum, as ASVs from Phycisphaeraceae did not significantly decrease in abundance in the other two Symbiodinium species. Kellogg (2019) found three OTUs belonging to Phycisphaerales to be bacterial associates across seven deep-sea corals (Eloe et al., 2011; Allers et al., 2013). Sequences similar to those three OTUs were not associated with tropical corals, possibly suggesting a role specific to cold-water corals (Kellogg, 2019). Contrarily, members of Phycisphaerales were associated with tropical stony corals, and in greater abundance under highly variable environmental conditions by Ziegler et al. (2017). Together with our results, these findings suggest that the observed decrease in Phycisphaeraceae with temperature may not be indicative of thermotolerance, but rather a consequence of the increase in relative abundance of Labrenzia in our most thermotolerant host.

In total, no ASVs were indicative of thermotolerance across all Symbiodinium species. Rather, we found ASVs that were differentially abundant by temperature differed by Symbiodinium species. These differences may be explained by functional redundancy of thermotolerance within and between bacterial communities (Stilianos et al., 2018). While Symbiodinium-associated bacterial communities may differ substantially, even among similar Symbiodinium species, these communities are likely to converge on key functional traits. For example, members of Rhodobacteraceae have generally been found to associate with other algae and provide vitamins and other amino acids to their hosts (Rosales et al., 2020), and play important roles in marine sulfur and nitrogen cycles (Chen et al., 2011). In microbial communities, many species fulfilling similar functions often coexist (Louca et al., 2018). Since all our cultures contained diverse members of Rhodobacteraceae, functional stability with stress could theoretically be maintained despite each culture containing different representatives of this family. For example, soil microbial communities showed high functional stability in long-term drought (Yuste et al., 2014). Future research on the relationship between phylogenetic diversity and functional stability may evidence similar consistencies in functional stability of dinoflagellate microbial communities under long-term increases in heat stress driven by climate change.

Additionally, high rates of horizontal gene transfer in long-term co-cultures could result in changes in the functional roles of different taxa, including those that may affect Symbiodinium thermotolerance (Cairns et al., 2018). This could have led to the development of similar thermotolerance-related functions in different taxa, leading to dissimilar patterns in bacterial community change with heat stress across cultures. While shifts in the functional roles of different community members have not been well-documented in long-term co-cultures, research on the human gut microbiome suggests that high rates of horizontal gene transfer and selective pressure have promoted functional redundancy in gut microbial communities despite wide diversity in microbiome community composition (Tian et al., 2020). More similar to our study system, biofilm bacterial communities of the kelp Ecklonia radiata commonly had horizontal gene transfers between bacterial members of the same class or order, and these transfers often included genes involved in stress response (Song et al., 2021). More research on Symbiodiniaceae is necessary to determine whether functional redundancy may be important for thermotolerance rather than any one bacterial taxa in particular.