Characterization of macroalgal-associated microbial communities from shallow to mesophotic depths at Manawai, Papahānaumokuākea Marine National Monument, Hawai‘i

Sampling site and abiotic factors

Samples were collected from August 03–08, 2019 from nine sites around Manawai PMNM, NWHI (Fig. 1). Field collections were approved by the Papahānaumokuākea Marine National Monument, permit number PMNM-2018-029. Collection numbers varied per site ranging from 6 to 12 samples (n = 6 at sites A and G; n = 8 at site H; n = 9 at sites B through F; n = 12 at site I). All sites varied by depth, except sites C and E. We categorized sampling into three depth zones according to Hinderstein et al. (2010): shallow subtidal (1.5, 2.0 m), subtidal (13, 22.5, and 27 m), and upper mesophotic (55, 58, and 75 m).

Light attenuation and temperature data for Manawai were characterized for each depth. Temperature data were collected using a Shearwater Petrel (Richmond, BC, Canada) wrist-mounted dive computer (upper and lower mesophotic) and Aqualung i200 (Vista, CA, USA) wrist-mounted diver computer (shallow subtidal and subtidal), which were calibrated prior to sample collection. Several minutes of immersion at the target depths were allowed for the dive computers to equilibrate with ambient temperature.

Six irradiance profiles were collected from the forereef and backreef at Manawai from July 15–20, 2021. The irradiance profiles were calculated as previously described (Spalding, Foster & Heine, 2003). Underwater irradiance was measured by lowering a spherical (4π) quantum sensor (Underwater LI-193SA; LI-COR, Lincoln, NE, USA) through the water column and storing the data in a datalogger (model 1400; LI-COR, Lincoln, NE, USA). The sensor was mounted on a lowering frame (LI-COR, Lincoln, NE, USA) and attached to the datalogger with a 20-m cable marked at 0.5-m intervals. PAR (μmol photons·m⁻²·s⁻¹) was recorded every 0.5 to 1 m to a maximum depth of 15 m. Profiles were completed from the sunny side of the boat on clear days between 1100 and 1300 h. K₀ was calculated from irradiance profiles according to Beer’s Law. The %SI (surface irradiance) was calculated from irradiance extrapolated from K₀ at 0.01 m.

Sample collection

Approximately one to three cm of new growth (i.e., the apical ends) of each most abundant macroalgal species (n = 3 per morphospecies) per site was collected. The most abundant macroalgal species were visually determined by experienced divers and phycologists with expertise on Hawaiian macroalgae. Each individual sample was collected into individual Whirlpak^© bags by the scientific dive team and sealed after collection underwater to limit any possible contamination among samples. Each sample was rinsed with 3.5% sterile artificial seawater to remove loosely attached epibionts and sand as previously described (Kuba et al., 2021). Individual thalli were then placed in RNAlater and stored overnight at 4 °C before freezing at −20 °C on the NOAA ship Rainier. A background seawater sample (n = 1) for each site was collected by filtering 50 mL of ambient seawater from each site (Boström et al., 2004) with an additional 50 mL artificial seawater control (n = 1) through a sterile 0.2 μm filter and preserved in 5 mL of RNAlater. This artificial seawater control was the same water used to rinse the samples prior to placing in RNAlater and served as a control for decontamination processing after assigning taxonomy. Each individual water sample was also collected into individual containers prior to filtering. The background seawater samples were collected approximately 1 m above the bottom at each site, and care was taken to collect water up current with no mixing with the benthos and represents a seawater control. These water controls serve as background controls to better understand the macroalgal-specific associations at our sample locations. This volume of seawater has an associated average DNA extraction efficiency of 92% for marine planktobacteria (Boström et al., 2004). Filters were stored at 4 °C before freezing at −20 °C. All samples for DNA analysis were shipped overnight on Techni Ice™ frozen at −80 °C from Honolulu, HI to Charleston, SC and subsequently frozen at −80 °C.

Identification of macroalgae

Macroalgae were identified visually based on morphological descriptions using Abbott (1999), Abbott & Huisman (2004), and Huisman, Abbott & Smith (2007) or by molecular analysis (Fig. 2; Table 1). Molecular analyses were completed and compared to voucher identifications (Table S1). For macroalgae that could not be conclusively identified by morphology, genomic DNA was extracted using an OMEGA E.Z.N.A.^® Plant DNA DS Kit (OMEGA Biotek, Norcross, GA, USA) or a NucleoSpin Plant II Kit (Macherey-Nagel, Düren, Germany) following the manufacturer’s protocol. For identification of the Chlorophyta specimens (vouchers NWHI 878, NWHI 1071 and NWHI 1049), a portion of the tufA gene (elongation factor Tu) was amplified using primers tufA_alg_up and tufA_alg_do (Handeler et al., 2010). For the Rhodophyta (NWHI 1073, NWHI 1032 and NWHI 813), a portion from the DNA barcode region near the 5′ end of the mitochondrial COI (cytochrome c oxidase subunit 1) gene was generated using primers GazF1 and GazR1 (Saunders, 2005). For further identification of the Rhodophyta (NWHI 1073, NWHI 1032 and NWHI 813), the rbcL (ribulose-1, 5-bisphosphate carboxylase/oxygenase large subunit) gene was amplified as two overlapping fragments using the primer pairs rbcLF7 and R898 (Gavio & Fredericq, 2002; Kim, Kim & Nelson, 2010) and rbcLF762 and R1381 (Kim, Kim & Nelson, 2010). Amplification conditions consisted of 94 °C for 3 min followed by 35 cycles of 30 s at 94 °C, 30 s annealing at 55 °C and 5 min synthesis at 72 °C, followed by a final extension at 72 °C for 7 min for tufA, and 4 min at 96 °C for denaturation, followed by 35 cycles of 60 s at 94 °C, 60 s at 42 °C and 90 s at 72 °C, with a final 10 min extension cycle at 72 °C and soak cycle at 10 °C for rbcL. The barcode region was amplified as previously described (Saunders, 2005). Successful PCR products were cleaned with ExoSAP-IT™ Express PCR Product Cleanup kit and submitted for sequencing at GENEWIZ (South Plainfield, NJ, USA). Raw sequence reads for each gene were assembled, edited, and aligned using the MUSCLE v. 3.8.425 plug-in (Edgar, 2004) in Geneious Prime 2021.0.3 (Fukunaga, 2008). Molecular sequences generated for the six macroalgae species were compared to those in GenBank using BLAST (Basic Local Alignment Tool; www.ncbi.nlm.nih.gov).

DNA extraction, library preparation, sequencing

All samples were thawed on ice and DNA was extracted from individual thalli using the FastDNA SPIN Kit (MP Biomedicals, Santa Ana, CA, USA) following the manufacturer’s protocol with minor modifications. Approximately 0.5 g of each algal sample was weighed into the lysis matrix E tube using ethanol flame-sterilized forceps. Thalli were split into two lysis tubes if an individual was >0.5 g. Each seawater control filter was cut and divided into lysis matrix E tubes to a final weight of 0.5 g. Extractions were completed on entire algal samples to include both associated epibionts and endobionts. Lysis tubes were placed in a cold aluminum rack and homogenized at 3,800 RPM for 30 s (BioSpec BeadBeater, Bartlesville, OK, USA). Bead beating was repeated twice with an incubation period of 30 s on ice in between each homogenization. DNA was eluted twice with 50 µl 0.1 mM Tris (pH 8.0). If thalli were split, eluted DNA was combined into one tube before further processing. DNA was quantified with a Qubit 3.0 fluorometer using the dsDNA high sensitivity kit (ThermoFisher Scientific, Waltham, MA, USA).

PCR was performed using the protocol from Klindworth et al. (2013) targeting variable regions 3 and 4 of the bacterial SSU rRNA gene as previously described (Kuba et al., 2021). Negative controls were included through PCR protocols, however no blank extractions were included. Primers possessed Illumina overhang sequences that were used for ligation of index sequences for all macroalgal and seawater controls. Cleaned amplicons were indexed and sequenced on an Illumina MiSeq per the manufacturer’s protocol, generating 2X300 base pair (bp) paired-end reads.

Bioinformatics of macroalgae-associated bacterial communities

Demultiplexed sequences with adapters removed were analyzed as previously described (Kuba et al., 2021). Amplicon sequence variants (ASVs) were inferred per MiSeq run (Callahan et al., 2016). Chimera identification and removal was performed after runs were merged using the “removeBimeraDenovo” command as implemented in dada2 (version 1.18). Each biological triplicate was kept as an independent sample for further processing to compare variation among individuals.

Statistical analyses were performed using R version 4.0.1 (R Core Team, 2010) and visualizations were generated using ggplot2 (Wickham, 2016), phyloseq (McMurdie & Holmes, 2013), and microbiome (Lahti & Shetty, 2017). Differential abundance analyses were utilized using the DESeq2 tool (Love, Huber & Anders, 2014) with an alpha <0.01 to identify the amplicon sequence variants (ASVs) contributing to overall differences among samples. This tool accounts for low dispersion estimates and is consistent across studies with varying replicates (Love, Huber & Anders, 2014). Data were then normalized using variance stabilizing transformation implemented with DESeq2 to compare microbiota across samples. Rarefaction of the number of reads per sample were visualized using the “rarecurve” command in vegan version 2.0–4 (Oksanen et al., 2016). Read abundances were transformed prior to calculating Euclidean distances. Visualizations of hierarchical clustering were performed with variance stabilized Euclidean distances. This transformation was applied due to the low associated false positive rate (McMurdie & Holmes, 2014). Alpha diversity was estimated using multiple indices (Table S3), but we chose to report the Simpson’s diversity index because of its emphasis on species evenness rather than richness as compared to the Shannon-Weaver index and Faith’s phylogenetic diversity; Simpson’s diversity index values closer to zero have lower microbial diversity. An ANOVA was run on each individual diversity model and a post-hoc Kruskal-Wallis test was ran on each individual model to compare between algal phylum/seawater. Beta diversity was visualized using a non-multidimensional scaling (NMDS) plot based on the Bray-Curtis dissimilarity metric. Bray-Curtis dissimilarity was applied rather than VST-stabilized Euclidean distances to account for distance and microbial species identity. Statistical tests, ANOVA, PERMANOVA, PERMDISP2, and ANOSIM were performed in vegan (Oksanen et al., 2016). If groups were non-homogenous, an ANOSIM was performed where the R-value ranges from 0 to 1 representative of indistinguishable to well-separated communities, respectively (Clarke & Gorley, 2001). Taxonomic distributions of macroalgal-associated microbiota were visualized using a heat map through ggplot2 (Wickham, 2016) with the associated cutoff for ASV relative abundance set to >0.05% and presence in >3 samples.

Data availability

Sequence data are available through NCBI Sequence Read Archive (SRA) BioProject number PRJNA833318. Macroalgal sequences used for taxonomic identification are available with accession numbers OR066430–OR066437, OK448437, OK448460 (Table S1). Relevant code used for the bioinformatics work is available on GitHub (https://doi.org/10.5281/zenodo.7975190).

Sample collection and abiotic factors

The most abundant macroalgae at Manawai were collected at each site (Fig. 1; Table 1). Key characteristics of these species are noted in Table 1. The calcification levels described include uncalcified, lightly calcified, and calcified. The shallow subtidal (1.5 to 2 m) temperatures ranged from 26–27 °C, subtidal (13 to 27 m) temperatures from 25–27 °C, and upper mesophotic (55 to 75 m) from 22–27 °C. The average diffuse attenuation coefficient (K_o) was −0.118 m⁻¹ (±0.01 m⁻¹ SE). Shallow subtidal irradiance ranged from 1,003 to 1,064 μmol photons m⁻² s⁻¹, with 79% to 84% surface irradiance (SI). Subtidal irradiance ranged from 53 to 274 μmol photons m⁻² s⁻¹, with 4% to 22% SI. Lastly, the upper mesophotic was characterized by irradiance ranging from 0.18 to 1.93 μmol photons m⁻² s⁻¹ and 0.014% to 0.152% of SI. The average 10% optical depth was 19 m (2.3/K_o), and the 1% optical depth was 39 m (4.6/K_o) (Kirk, 1994). Light attenuation and temperature data are provided in Table S2.

Molecular identification of macroalgae

Four of the six macroalgal species that could not be identified based on morphology alone were identified to at least genus level using BLAST. The Chlorophyta macroalgal voucher NWHI 1049 was 100% identical to Umbraulva kaloakulau H.L. Spalding & A.R. Sherwood (Chlorophyta) sequences, while the other two Chlorophyta algal vouchers were both identified as Halimeda velasquezii W.R. Taylor (Chlorophyta), having 99.2% and 99.6% (NWHI-878 and NWHI-1071, respectively) similarity with H. velasquezii sequences in GenBank, as well as additional morphological confirmation. For molecular identification of the Rhodophyta, two markers were used: COI and rbcL, and combined with morphological analysis. BLAST results for voucher NWHI-813 indicated it was Laurencia galtsoffii M. Howe with 99.8% of similarity for the COI marker (amplification of rbcL was not successful for this specimen). COI sequence from voucher NWHI-1032 yielded 93.2% similarity with Gracilaria shimodensis (Rhodophyta) and was 96.3% to G. hayi, while COI sequence from voucher NWHI-1073 was only 85.8% similar to Symphyocladiella gracilis (Rhodophyta) and rbcL sequence was 93.5% similar to Wrightiella tumanowiczii (Rhodophyta). Because additional morphological analysis is necessary to attribute species names to these vouchers, we referred to these specimens as Gracilaria sp. and Wrightiella sp.

Characterization of macroalgal associated microbial communities at Manawai

A total of 77 macroalgal specimens representing 15 species, nine background seawater samples, and one artificial seawater control were sequenced, resulting in 5,357,308 total reads representing 50,445 ASVs with an average length of 415 ± 30 bp after quality control, filtering, and chimera filtering. Sequencing appeared to capture the overall diversity for each specimen based on the rarefaction curve plateau across all samples (Fig. S1). Sample read depth varied while characterizing microbial diversity (medium, 18,722; median, 165,447; maximum, 364,285 across samples).

Overall the average microbial richness was highest for the Chlorophyta samples followed by Rhodophyta and then Ochrophyta, as measured by number of observed ASVs, Shannon index, Simpson’s index, and Faith’s phylogenetic diversity (Fig. 3 and Fig. S2). The number of observed ASVs for the background seawater was lower than all macroalgal samples (Fig. 3). Microbial communities associated with macroalgae at Manawai significantly overlapped (ANOSIM R = 0.1955; p = 0.0039) (Fig. S3). The associated phylogenetic diversity was significantly different between Ochrophyta and both Chlorophyta and Rhodophyta (Fig. 3 and Fig. S2; Tables S3 and S4). According to the Simpson’s index, microbial diversity was lowest for Ochrophyta (Fig. 3). The total number of observed ASVs across all Ochrophyta specimens ranged from 229 (Distromium sp.) to 5,942 (Dictyota ceylanica Kützing) both from site E at 55 m (Table S3). The number of observed ASVs varied among species and among biological triplicates. A specimen of H. velasquezii (Chlorophyta) collected at 2 m depth had the lowest diversity based on the Simpson’s index (0.078), followed by another replicate of the same species at the same depth and site (0.106).

Actinobacteriota, Alphaproteobacteria, Bacteroidota, Cyanophyceae, Gammaproteobacteria, Planctomycetota, and Verrucomicrobiota were the most relatively abundant bacterial classes associated with all macroalgal phyla (Fig. 4). Myxococcota was not associated with Ochrophyta, and Firmicutes was unique to Chlorophyta. The background seawater samples varied in associated bacterial taxa, but Gammaproteobacteria sequences had the highest relative abundance across all controls (Fig. 4). Ochrophyta were associated with fewer bacterial taxa as compared to the other phyla and the seawater controls. Gammaproteobacteria associated with Ochrophyta (Pseudomonadales and Alteromonadales) showed the highest relative abundance.

Hierarchical clustering based on Euclidean distances separated the microbial communities into three groups (Fig. 5). Each of these clusters contained multiple phyla and species; however, only Rhodophyta and Chlorophyta were represented in cluster one (Fig. 5). Throughout the three clusters, communities from each macroalgal species clustered together, with only a few exceptions (Fig. 5). U. kaloakulau from 75 m depth was found in two different clusters, with similar microbial communities to Ochrophyta and Rhodophyta species. The seawater controls grouped together in cluster two regardless of collection site (Fig. 5). Microbial community structure from the same macroalgal species was similar. All Ochrophyta species and specimens were clustered in the second cluster.

Spatial comparison of M. setchellianum microbial community structure

M. setchellianum was the most abundant macroalga at Manawai, which resulted in its collection from six sites spanning all depth zones. This replication allowed for the examination of microbial communities over depths and depth zones of a single macroalgal species. Based upon overall variation across sites, samples were analyzed by depth zone rather than absolute depth. Specimens of M. setchellianum collected at sites A (1.5 m) and I (27 m) were associated with the highest diversity (>0.99). Most M. setchellianum communities were found within cluster three of the hierarchical clustering dendrogram (Fig. 5). Samples from sites C and F (n = 3 at both sites), both from the upper mesophotic zone, were found in a subgroup of cluster three that only contain macroalgae sampled from the upper mesophotic (Fig. 5). All seawater were found in a subgroup of cluster two rather than as an outgroup, indicating some similarity between microbiota associated with macroalgal hosts and surrounding seawater (Fig. 5). The M. setchellianum communities from sites A, B, and H (n = 3 at all sites) from the shallow subtidal and subtidal zones form the other subgroup of cluster three, along with Wrightiella sp. (also from the subtidal). The outliers were the subtidal M. setchellianum communities from site I, which were found within cluster one (Fig. 5). Microbial communities associated with M. setchellianum overlapped significantly across sites (ANOSIM R = 0.38; p = 0.0049) and based on collection depth zone (ANOSIM R = 0.36; p = 0.0059). Shallow subtidal M. setchellianum microbial communities appeared to differ the most as confidence intervals for subtidal and upper mesophotic communities overlapped (ANOSIM R = 0.11; p = 0.002) (Fig. 6A).

Microbial community structure of Halimeda spp.

Two Halimeda species (H. discoidea Decaisne and H. velasquezii) were abundant at three shallow subtidal and subtidal sites and co-occurred at one site (I, 22.5 m depth). All Halimeda spp. bacterial communities were within cluster one except for H. velasquezii bacteria from site D from the shallow subtidal (Fig. 5). Microbial associations were significantly different based on site when each species was examined individually (PERMANOVA p = 0.003 for H. discoidea, p = 0.007 for H. velasquezii) (Fig. 6B). Close associations were also apparent among the biological triplicates for both Halimeda spp. at each site (Figs. 5, 6B). Analysis by PERMANOVA supported a difference in overall microbial communities associated with Halimeda spp. based on site (p = 0.001) (Fig. 5B). Bacterial communities associated with H. discoidea and H. velasquezii at all sites were different but overlapped (ANOSIM R = 0.56; p = 0.0001) (Fig. 6B). Microbial communities were well separated between these species at the single site where they co-occurred, however, this separation was not significant (ANOSIM R = 0.82, p = 0.1) (Fig. 6C).

Microbial communities associated with the cryptogenic alga C. tumulosa

C. tumulosa was collected from two sites (B and D, n = 3 at each site). At site D (2 m), we also collected specimens of native rhodophyte L. galtsoffii (n = 3), allowing us to examine differences in microbial communities between a cryptogenic and native alga in the same family (Rhodomelaceae). Specimens of C. tumulosa and L. galtsoffii were associated with five overlapping dominant associated bacterial class, similar to other rhodophycean algae sampled (Fig. S2B). However, C. tumulosa was associated with two unique dominant bacterial classes: Gammaproteobacteria and Cyanobacteriia (Fig. 4B). The Simpson’s diversity index associated with all C. tumulosa averaged 0.959 and L. galtsoffii averaged 0.963 (Table S3). All bacterial communities from C. tumulosa specimens were in cluster one in the hierarchical clustering analysis but in separate subgroups, independent of sampling site (Fig. 5). Communities from L. galtsoffii were also found in cluster one, but in a different subgroup (Fig. 5). Microbial community composition from these two host species were separated without statistical significance (ANOSIM R = 1; p = 0.1, Fig. S4).