Colonies of the soft coral E. rubra were collected at a depth of approximately 15 to 20 m near Eoyudo, located off the southern coast of Korea (34°39′14.00′′N, 128°34′19.00′′E), in November, 2023 using standard SCUBA techniques, and the seawater temperature during collection was 16–17°C ( Figure 1 ). The specimens were transported to a laboratory aquaria facility and acclimated for three days in filtered seawater maintained at 16°C and 33 ppt salinity, under a 14-hour light and 10-hour dark cycle at the ambient light intensity (800–900 lux). Seawater, pumped from a depth of 10 m in proximity to the institute (300 m from the coast) and filtered serially through 100, 10, and 1 µm filters (Filtertech, South Korea), was continuously supplied to four 100 L tanks at a minimal rate of 50 L h⁻¹ for the coral acclimation, and seawater finally filtered through a 0.1 µm membrane was used for the thermal stress experiment. After acclimation, three intact colonies of corals were assigned to 26°C, 28°C, and 30°C experimental tanks respectively for 24 hours. The control group for the thermal experiment was kept in seawater at 16°C, and the experimental conditions were controlled by AquaController III Pro (Neptune Systems, Morgan Hill, CA, USA). Upon completion of the exposure, coral samples were immediately stored in an ultra-low temperature freezer at –80°C until DNA extraction.

Frozen tissue samples were ground in a mortar with liquid nitrogen. Approximately 200 mg of powdered tissue was homogenized in 1 mL of CTAB lysis buffer containing 3% CTAB, 5.0 M NaCl, 3% PVP, 0.5 M EDTA (pH 8.0), 1 M Tris-HCl (pH 8.0), and 0.2% β-mercaptoethanol. The homogenate was incubated at 65°C for 1 hour. An equal volume of phenol:chloroform:isoamyl alcohol (25:24:1) was then added to denature proteins, followed by centrifugation at 10,000 ×g for 15 minutes at room temperature. The aqueous phase was collected, and RNase A (30 mg/mL) was added, followed by a 1-hour incubation at 37°C. DNA was further extracted with phenol:chloroform:isoamyl alcohol (25:24:1), then purified again using chloroform:isoamyl alcohol (24:1). After centrifugation at 10,000 ×g for 15 minutes at room temperature, the aqueous phase was recovered. DNA was precipitated by adding one-tenth volume of 3 M sodium acetate (pH 5.2) and an equal volume of 100% ethanol. The precipitated DNA was washed with 70% ethanol and resuspended in ion-exchanged ultrapure water. DNA concentration was measured using the Picogreen method with a Victor 3 fluorometer (Perkin Elmer Inc., Waltham, MA, USA), and quality was confirmed by agarose gel electrophoresis.

Sequencing libraries were prepared according to the Illumina 16S Metagenomic Sequencing Library protocol to amplify the V3–V4 region. A total of 2 ng of genomic DNA was amplified in a PCR reaction containing 5× reaction buffer, 1 mM dNTP mix, 500 nM of each universal forward and reverse primer, and Herculase II fusion DNA polymerase (Agilent Technologies, Santa Clara, CA, USA). The first PCR was conducted with an initial denaturation at 95°C for 3 minutes, followed by 25 cycles of 30 seconds at 95°C, 30 seconds at 55°C, and 30 seconds at 72°C, and a final extension at 72°C for 5 minutes. The primers used in the first amplification included Illumina adapter overhangs: V3-F (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3^′) and V4-R (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3′). The resulting PCR product was purified using AMPure beads (Agencourt Bioscience, Beverly, MA, USA), and 2 µL of the purified product was used as the template for a second PCR using NexteraXT indexed primers. The second PCR followed the same thermal cycling conditions as the first but was run for 10 cycles. The final PCR products were again purified using AMPure beads. Quantification was performed using qPCR (KAPA Library Quantification Kit for Illumina platforms), and quality was assessed using the TapeStation D1000 ScreenTape system (Agilent Technologies, Waldbronn, Germany). Paired-end sequencing (2 × 300 bp) was conducted on the Illumina MiSeq platform by Macrogen (Daejeon, South Korea).

After sequencing, the raw MiSeq data were demultiplexed based on index sequences, and fastq files were generated for each sample (Chen et al., 2018). Adapter sequences were removed using the fastp program, and sequencing errors in the overlapping region of paired-end reads were corrected. Paired-end reads were merged into single sequences using FLASH (v1.2.11) (Magoč and Salzberg, 2011). Assembled reads shorter than 400 bp or longer than 500 bp were removed. Low-quality sequences, ambiguous bases, and chimeric reads were filtered out using CD-HIT-OUT (Fu et al., 2012), and the remaining reads were clustered into species-level operational taxonomic units (OTUs). Representative sequences for each OTU were identified using BLAST+ (v2.9.0) (Zhang et al., 2000) against the NCBI 16S Microbial Reference Database. Taxonomic assignments were based on the best hit, but OTUs were left unclassified if the query coverage or identity of the best hit was below 85%. The resulting OTU table was used for comparative analyses of microbial communities using QIIME (v1.9) (Caporaso et al., 2010).

Microbial diversity within samples was evaluated using alpha diversity indices, including Shannon, Inverse Simpson, and Chao1, calculated based on rarefied OTU tables. Comparisons between the control and heat stress groups were performed using Student’s t-test and visualized as boxplots with Plotly. To assess beta diversity, principal coordinate analysis (PCoA) was performed using both weighted and unweighted UniFrac distances. To ensure direct comparability with the functional data, we also calculated Bray–Curtis dissimilarities on the rarefied OTU relative abundance table (square root transformed). The resulting distance matrix was subjected to principal coordinates analysis (PCoA) and to PERMANOVA and PERMDISP with 999 permutations using scikit-bio ( Supplementary Table S4 ). Group distribution trends were visualized with 95% confidence ellipses using the EllipticEnvelope method. To verify sequencing depth sufficiency, rarefaction curves were generated using Chao1 and Observed OTUs. All samples showed curves approaching saturation, indicating adequate coverage for diversity estimates ( Supplementary Figure S1 ).

Taxonomic composition was compared at the genus and phylum levels based on relative abundance and visualized as stacked bar plots. Genus-level composition was further evaluated using rank–abundance curves to simultaneously assess dominant and rare taxa. To statistically compare phylum-level relative abundances between control and 26°C-treated (heat stress) groups, both Welch’s t-test and the Mann–Whitney U test were applied to account for variance and non-normal distribution. The average difference (Δ) was calculated as the mean relative abundance difference between groups ( Supplementary Tables S1 , S2 ). Functional analysis was performed using PICRUSt2 to predict KEGG and MetaCyc pathways. Significant changes were evaluated based on log₂ fold change and the Mann–Whitney U test. Pathways with marked increases or decreases were visualized using bar plots, and the top 50 most abundant functions were shown as a log-scaled heatmap. Functional PCoA based on Bray–Curtis distances was also conducted and visualized with Seaborn and Matplotlib. LEfSe analysis was used to identify significantly enriched species-level taxa, and the results were visualized using a bidirectional bar plot based on log-transformed LDA scores. Lastly, to evaluate correlations between genus-level taxa and predicted metabolic functions, Spearman correlation analysis was performed, and relationships among the top 20 genera and functions were visualized as a heatmap. All analyses were conducted using various Python libraries, including Pandas, Scipy, Seaborn,Matplotlib, andPlotly (Barrett et al., 2005; McKinney, 2011; Sievert, 2020; Virtanen et al., 2020; Waskom, 2021).

To investigate changes in the symbiotic microbial community of corals under elevated seawater temperatures, colonies of E. rubra (a soft coral species) were collected from coastal waters near Eoyudo, located off the southern coast of Korea. The collected corals were transferred to experimental aquaria and acclimated for three days in filtered seawater maintained at 16°C. Based on the sea surface temperature data around South Sea, Korea (Ribas-Deulofeu et al., 2023) and the other soft coral, S. gracillimum, study on heat stress in our laboratory (Woo and Yum, 2022), the heat stress experiment of E. rubra was designed originally at 26°C, 28°C, and 30°C but the groups of 28°C and 30°C were melted down within 24 hr and only 26°C group was compared with the control group to investigate its microbial community. Total DNA was extracted from each sample, followed by 16S rRNA-based metagenomic sequencing.

The impact of thermal exposure on the microbial community of E. rubra was assessed through alpha and beta diversity analyses. All alpha diversity indices (Shannon, Chao1, and Inverse Simpson) showed a modest increase under the heat stress condition, indicating a potential expansion of microbial diversity in response to heat stress ( Figure 2 ). The similar alpha diversity values suggest that the overall microbial diversity remained relatively stable. However, high inter-individual variability precluded statistical significance. To evaluate whether this non-significance was simply due to small sample size (n = 3 per group), we performed a post-hoc power analysis using Cohen’s d for each index. Moderate-to-large effect sizes were observed for Shannon (d =0.79) and Inverse Simpson (d =0.93), yet the resulting statistical power was low across all metrics (Shannon ≈11.7%, Inverse Simpson ≈14.4%, Chao1 ≈5.2%; see Supplementary Table S3 ).

For beta diversity, principal coordinate analysis (PCoA) was performed using both weighted and unweighted UniFrac distances ( Figures 3a, b ). In the weighted UniFrac analysis, PC1 explained 72.8% of the total variance, and the centroids of the two groups were located in close proximity, indicating that the composition and relative proportions of dominant taxa remained largely unchanged after heat stress. In contrast, the unweighted UniFrac analysis showed greater separation between groups, with PC1 and PC2 explaining 25.1% and 21.9% of the variance, respectively. In addition to the UniFrac-based comparisons, we conducted a complementary Bray–Curtis β-diversity analysis on the same OTU relative-abundance table. Using Bray–Curtis distances, the first two PCoA axes accounted for 45.1% and 21.8% of the variance, respectively ( Figure 3 ). Although PERMANOVA (pseudo-F =1.08, p =0.205) and PERMDISP (F =0.099, p =0.786) did not indicate statistically significant group separation ( Supplementary Table S4 ), the ordination nonetheless displayed a visibly clearer segregation between CONTROL and HEAT-STRESS samples than the UniFrac plots.

This pattern was also reflected in the rank-abundance curve analysis ( Figure 4b ). Compared to the control group, the heat stress group exhibited higher overall abundance, and the curve was more gradual in the latter half (after rank 150), indicating increased evenness and a relative enrichment of rare taxa. These results are consistent with the alpha and beta diversity findings, supporting the interpretation that heat stress may induce community restructuring through increasing rare microbial populations.

To evaluate the effects of heat stress on the symbiotic microbial community structure of E. rubra in more detail, we analyzed the distribution and relative abundance of bacterial taxa. Among a total of 599 species identified across both groups, 389 species (64.9%) were shared between the control and heat stress samples, while 77 species (12.9%) were unique to the control group and 133 species (22.2%) were detected only under heat stress conditions ( Figure 4a ). Taxonomic composition was further examined at both the phylum and genus levels ( Figures 4c, d ; Supplementary Tables S1 , S2 ). At the phylum level, the dominant group Proteobacteria increased from an average of 41.7 ± 8.6% in the control group to 50.3 ± 7.3% in the heat stress group. Bacteroidetes also showed an increasing trend, from 24.0 ± 5.7% to 29.1 ± 7.5%, although these changes were not statistically significant (p > 0.05). Several minor phyla, including Acidobacteria, Firmicutes, and Nitrospinae, also showed increasing trends under heat stress, but without statistical significance. At the genus level, Alteromonas, Phaeobacter, Pseudomaricurvus, Oleibacter, and Pseudoruegeria were significantly enriched in the heat-stress group (p < 0.05), while Thiohalophilus, along with flavobacterial genera such as Zobellia and Paracoccus, showed significant decreases in relative abundance (p < 0.05).

Given the high variation among samples, group-level comparisons based solely on average relative abundance can be misleading. To better quantify taxonomic shifts associated with heat stress, LEfSe (Linear Discriminant Analysis Effect Size) analysis was conducted ( Figure 5 ; Supplementary Table S5 ). Using a threshold of LDA score > 2.0 and p < 0.05, a total of 43 species showed significant differences between the groups. The top five species enriched under heat stress were Colwellia aquamaris, Wenzhouxiangella marina, Microbulbifer hydrolyticus, Poseidonibacter lekithochrous, and Pelagimonas phthalicica. In contrast, species that were less abundant in the heat stress group included Endozoicomonas elysicola, Thiohalocapsa marina, and Dulcicalothrix necridiiformans, some of which are known for their photosynthetic capacity or benthic adaptations. For completeness, we also repeated the LEfSe analysis using a stricter LDA cut-off of 3.0. As a result, only 15 discriminatory species were retained, all of which are a subset of those shown in Figure 5 . This result is provided in Supplementary Figure S2 .

Changes at the phylogenetic level alone do not provide direct insight into how microbial community restructuring may influence the physiological functions of E. rubra. In marine microbial ecosystems, taxonomically distinct organisms can share similar metabolic pathways due to functional redundancy. Therefore, an increase or decrease in the abundance of a specific taxon does not necessarily translate into changes in overall metabolic capacity. Conversely, even minor changes in taxonomic composition may lead to substantial shifts in host–microbe interactions if the expression of key metabolic pathways is altered. To address this, we reconstructed functional profiles based on predicted metabolic potential using PICRUSt2 and examined how temperature changes affected the overall metabolic spectrum of the community. We also investigated how these functional shifts correlated with specific bacterial taxa. A total of 366 metabolic pathways were predicted and quantified using the MetaCyc database in PICRUSt2, and Bray–Curtis distance-based PCoA was conducted. PC1 explained 68% of the total variance, and samples from the control and heat stress groups were separated ( Figure 6a ). A heatmap based on the log₁₀-transformed median abundance of predicted metabolic pathways also revealed overall differences in functional profiles between the groups ( Figure 6b ), which is consistent with the separation pattern observed in the PCoA analysis. Based on the median abundance of each pathway, log₂ fold changes were calculated. Among the 366 predicted pathways, 47 (12.8%) showed at least a twofold increase (log₂FC > 1), and 25 (6.8%) exhibited at least a twofold decrease (log₂FC < –1) under heat stress ( Supplementary Table S6 ). The syringate degradation pathway showed the largest increase, followed by elevated levels of pathways related to carbohydrate and cofactor metabolism, including enterobacterial common antigen biosynthesis, methylaspartate cycle, adenosylcobalamin biosynthesis, and isoprene biosynthesis. In contrast, pathways involved in membrane components and fermentation, such as vitamin E biosynthesis, chondroitin sulfate degradation, and pyruvate fermentation to acetone, were reduced. These results suggest that metabolic pathways may be reorganized in response to thermal stress ( Figures 7a, b ). To evaluate how these functional changes were linked to taxonomic composition, we calculated Spearman correlation coefficients between the top 20 most abundant genera and the predicted metabolic functions, and visualized the results in a heatmap ( Figure 8 ). Genera such as Thalassotalea, Amphritea, and Rubritalea exhibited positive correlations with most major metabolic pathways, especially those involved in amino acid and nucleotide biosynthesis. In contrast, genera such as Alteromonas and Nitrospira showed negative correlations with the same pathways, indicating that different genera contributed differentially to the functional profiles under heat stress. Interestingly, we observed cases where changes in genus-level abundance did not align with changes in functional output. This suggests that the taxa responsible for particular metabolic functions may have shifted under heat stress, even when overall abundance patterns remained relatively stable.