Examining the giant barrel sponge species complex: molecular and microbial differentiation of Xestospongia testudinaria in Singapore

The giant barrel sponges (Xestospongia spp.) belong to a pan-global species complex with evidence suggesting they could encompass up to 9 cryptic species. In this study, we leveraged molecular and microbial techniques to investigate giant barrel sponges (X. testudinaria) from Singapore in relation to their placement within this species complex. Twenty-four giant barrel sponges from three sites were sequenced with mitochondrial (CO1) and nuclear (ATP6) DNA markers, identifying 6 distinct haplotypes belonging to 4 of the proposed barrel sponge species. Analysis of the X. testudinaria microbiomes was achieved with 16S rRNA gene amplicon sequencing. The microbiome composition of X. testudinaria did not differ by reef site, deviating from a pattern frequently observed in coral microbiomes across Singapore. However, there was significant differentiation in microbiome composition by host genetics consistent with the proposed species boundaries. General linear models identified 85 amplicon sequence variants (ASVs) as highly significant (P < 0.01) in differentiating among the four Species Groups, consisting of 12 Archaea and 73 Bacteria, with the largest representation from phylum Chloroflexi. We also identified 52 core ASVs present in all sponges representing 33.0% of the total sequence reads. Our results support previous findings of microbiome differentiation in co-occurring genetic haplotypes of barrel sponges from the Caribbean. Together these studies underline the potential for ecological partitioning based on genetic haplotype that could contribute to cryptic speciation within the giant barrel sponge species complex.


Introduction
Sponges represent a functionally important component of coral reef ecosystems (Pawlik and McMurray 2020).Through active filter feeding, sponges contribute to ecosystem carbon (Perea-Blázquez et al. 2012;de Goeij et al. 2013;Pawlik and McMurray 2020), nitrogen (Wilkinson and Fay 1979;Hoffman et al. 2009;Fiore et al. 2013a, b), phosphorus (Zhang et al. 2015), and sulfur cycling (Zhang et al. 2019).These pathways are mediated by the sponge microbial symbionts (Morganti et al. 2017;Zhang et al. 2019;Rix et al. 2020).Coral reef degradation and subsequent reductions in coral cover (Wolff et al. 2018;Eddy et al. 2021) result in a dynamic system of species interactions between, corals, algae, and sponges (González-Rivero et al. 2011;Pawlik and McMurray 2020).Many reefs are experiencing an increase in sponge abundance (McMurray et al. 2010;Bell et al. 2013;Powell et al. 2014;Helber et al. 2018;Marlow et al. 2019).Thus, the weighted ecological role of reef sponges in ecosystem functioning could have greater importance as their relative abundance in these systems increases.
The giant barrel sponges (Xestospongia spp.) are ubiquitous on tropical reefs in both the Atlantic and Indo-Pacific (McMurray et al. 2010;Swierts et al. 2017;McGrath et al. 2018).Barrel sponges can grow well over a meter, capable of processing large volumes of seawater (McMurray et al. 2008(McMurray et al. , 2014)), particularly via the uptake of dissolved organic carbon (DOC) and nitrogen cycling pathways (Fiore et al. 2013a;Fiore et al. 2015;McMurray et al. 2016).Originally thought to contain three species-X.muta (Schmidt 1870) in the Caribbean and X. testudinaria (Lamarck 1815) and X. berquistia (Fromont 1991) in the Indo-Pacific-the barrel 166 Page 2 of 8 sponges are now considered part of a pan-global species complex with 17 genetic clusters, potentially representing 9 cryptic species (Swierts et al. 2017).Cryptic speciation can mask potential differences in ecological functioning and niche separation (Xavier et al. 2010;Fišer et al. 2017).Therefore, it is important to examine the potential for cryptic speciation particularly in dominant and resilient reef megafauna, like the giant barrel sponges (Bell et al. 2014).
Sponges possess diverse and highly species-specific microbiomes (Taylor et al. 2007;Lee et al. 2011;Schmitt et al. 2012;Thomas et al. 2016); however, intraspecific differences can result from geographic separation (Fiore et al 2013b;Marino et al. 2017) and genetic differences (Griffiths et al. 2019;Díez-Vives et al. 2020;Easson et al. 2020).Indeed, within the Caribbean giant barrel sponges, microbiome differentiation has been detected among sympatric individuals (Evans et al. 2021) consistent with the species boundaries proposed by Swierts et al. (2017).Demographic analyses of the same population revealed a shift in the dominant genetic group (Deignan et al. 2018), which when coupled with the microbiome data, suggests that there may also be a shift in ecological functioning or niche partitioning occurring for these important reef inhabitants.In the Indo-Pacific, barrel sponge microbiomes were examined across regional scales and found to differentiate by geographic location (Swierts et al. 2018); however, no fine-scale analysis to link the microbiome with genetic haplotype has been conducted within a single location.In this study, we examined the relationship between genetic separation and microbiome composition of giant barrel sponges from within Singapore to investigate whether Indo-Pacific barrel sponge microbiomes differentiate following the proposed species boundaries, as observed in the Atlantic.

Haplotype assignment
Genomic DNA was extracted from the inner (endosome) sponge tissue of the cryopreserved X. testudinaria with QIAGEN's DNeasy PowerSoil extraction kit before using Zymo Research's Genomic DNA Clean & Concentrator kit to remove potential inhibitors.All 24 sponges were sequenced for mitochondrial cytochrome oxidase 1 (COI) and adenosine triphosphate synthase subunit 6 (ATP6) genes using primer set CI-J2165 (5′-GAA GTT TAT ATT TTA ATT TTACCDGG-3′), CI-Npor2760 (5′-TCT AGG TAA TCC AGC TAA ACC-3′) and ATP6porF (5′-GTA GTC CAG GAT AAT TTA GG-3′), ATP6porR (5′-GTT AAT AGA CAA AAT ACA TAA GCC TG-3′), respectively.The PCR protocol for ATP6 amplification followed an initial denaturation at 95 °C for 15 min and 35 cycles of denaturing at 94 °C for 30 s, annealing at 40 °C for 45 s and extension at 68 °C for 1.5 min, with a final extension of 72 °C for 10 min.COI amplification was done at an initial denaturation of 95 °C for 5 min and 35 cycles of denaturing at 94 °C for 30 s, annealing at 40 °C for 60 s and extension at 72 °C for 1.5 min, with a final extension of 72 °C for 10 min.PCR was conducted in triplicate with a total reaction volume of 25 μL containing 0.5 μL of each primer (10 μM), 1 μL BSA (Bovine Serum Albumin; 200 ng/μL), 10 μL HotStarTaq Plus Master Mix (Qiagen), 12 μL sterile water, and 1 μL DNA template (undiluted purified DNA of varying concentrations).Purity of the DNA was measured spectrophotometrically with NanoDrop™ 2000 (Thermo Scientific) and quantified with Qubit™ 2.0 fluorometer.PCR products were purified with Zymo Cleanup and Concentrator kit and sent to a commercial company for DNA sequencing via Applied Biosystems 3730xl DNA Analyzer after processing with BigDye Terminator v3.1 Cycle Sequencing Kit.
Sequence reads were edited, aligned, and concatenated using Geneious Prime ® 2020.2.The identity of Xestospongia sequences was confirmed via BLAST.Haplotype network was constructed with PopArt v1.7 (Leigh and Bryant 2015).The sequencing reads were uploaded to the NCBI Database for the ATP sequences (GenBank accession numbers: OM154214-37) and the COI sequences (GenBank accession numbers: OM154238-61).

Sponge microbiome analysis
DNA was extracted from both inner (endosome) and outer (ectosome) sponge tissues with QIAGEN's DNeasy Power-Soil extraction kit according to the manufacturer's protocols and purified with Zymo Research's Genomic DNA Clean & Concentrator kit.The extracted DNA was quality-checked using NanoDrop™ 2000 Spectrophotometer, and Qubit™ 2.0 Fluorometer with dsDNA broad-range assay kits.Differentiation of inner and outer tissues was determined by the tissue color-with a deeper reddish brown for the ectosome tissue samples and white for the endosome tissue samples.This resulted in 48 total samples from the 24 sponges.515F (Parada) and 806R (Apprill) primers were used to amplify the V4 region of the 16S rRNA gene (Caporaso et al. 2011;Apprill et al. 2015;Parada et al. 2016) for sequencing of the prokaryotic portion of the sponge microbiomes.PCR amplification was done in 20 μL volumes with the following recipe: 10 μL of HotStarTaq, 1 μL each of 515F and 806R primers (10 μM), 1 μL dimethyl sulfoxide, 1 μL BSA (Bovine Serum Albumin; 200 ng/μL), 4 μL sterile nucleasefree water, and 2 μL DNA template.Amplification was carried out using an initial denaturing step of 95 °C for 5 min followed by the following 35 cycles of: 94 °C for 30 s, 53 °C for 40 s, and 72 °C for 1 min; and a final elongation step of 72 °C for 10 min.Following PCR, the PureLink Quick Gel Extraction Kit (Invitrogen) was used to excise and purify the bands of interest.The purified DNA samples were then quality-checked on an Agilent TapeStation using Agilent D1000 ScreenTapes and reagents before being sent for amplicon sequencing on an Illumina MiSeq platform at the Singapore Centre for Environmental Life Sciences Engineering (SCELSE).The raw sequence reads were uploaded to the NCBI Sequence Read Archive under BioProject accession number PRJNA976080.
Amplicon sequence reads were processed using DADA2 version 1.16 to generate amplicon sequence variants (ASVs) for each sample.Sequence data were processed twice, once for the multi-tissue samples where both inner and outer tissue sample reads were included individually, and again for the combined tissue, wherein both inner and outer tissue sample reads were merged into a single read file for each sponge.Sequence reads were trimmed to 210 for the forward read and 165 for the reverse read before error learning algorithms were applied and chimeric sequences removed.Contaminating sequence reads based on comparison to blank extractions were removed using the decontam package (Davis et al. 2018).Taxonomy was assigned to the genus level using the built-in native Bayesian classifier assign taxonomy based on the SILVA SSU r138.1 database.Sequences identified as mitochondria, chloroplast, or unassigned to a Domain were removed.Finally, rarefaction curves were used to assess samples that reached their ASV maximum, and samples were rarefied to account for variation in sequencing depths for the calculation of diversity metrics (Weiss et al. 2017) to 97,655 sequence reads per sample and to 101,990 sequence reads per sponge for the merged analysis.
Permutational analyses of variance (PERMANOVA) based on square root transformed Bray-Curtis matrices were conducted in PRIMER v7 for comparisons between the following groups: inner and outer sponge tissue, site, haplotype, and Species Group.Follow-up permutational multivariate analyses of dispersion (PERMDISP) were performed on haplotype and Species Group comparisons to test for homogeneity of dispersion.Differences between groups were visualized with non-metric multi-dimensional scaling (nMDS) plots.Alpha diversity metrics were calculated using vegan and compared by Species Group using ANOVA.Core, defined as presence in all samples, and unique ASVs by Species Group were also determined.Finally, generalized linear models (GLM) with negative binomial distribution were performed using the mvabund package (Wang et al. 2012) to identify the top ASVs contributing to the differences among Species Groups.The dataset was subsampled to 500 ASVs for the GLM analysis to focus on the most abundant ASVs driving the group distinctions.

Results and discussion
24 COI and ATP6 sequences were combined and trimmed to a final length of 989 base pairs.A total of 6 haplotypes, encompassing 8 polymorphic sites, were found.All 6 haplotypes were previously identified by Swierts et al. (2017) within Singapore's St John's Island Xestospongia population, which found 7 haplotypes in a total of 15 samples.Therefore, the haplotypes are annotated following Swierts et al. (2017).C6A2 and C2A1 were the most prevalent haplotypes and were found across all sampling sites (Table 1).Two unique haplotypes, C5A4 and C1A1 were found at Hantu and Sister's respectively.However, both haplotypes were present at St John's Island during the 2017 study by Swierts et al.; therefore, their absence from other sites is likely an artifact of sample size and not genetic structure within the population.The 6 haplotypes identified encompassed 4 Species Groups, as proposed for X. testudinaria by Swierts et al. (2017;Table 1).
There was no significant difference in microbiome composition between the inner and outer tissue samples (Pseudo-F = 0.90733, P = 0.5196).Therefore, only the combined sequence data by individual X. testudinaria were used for further analysis.X. testudinaria microbiomes did not differ by site (Pseudo-F = 1.5765,P = 0.0543).Site differences in microbiome composition are common for corals within Singapore (Wainwright et al. 2019;Fong et al. 2020;Deignan and McDougald 2022;Moynihan et al. 2022;Deignan et al. 2023), highlighting how sponges relate to their microbiomes differently than other dominant benthic reef organisms.
Within the nMDS plot, the 6 haplotypes appeared to group together into 4 clusters, with haplotype C1A1 (n = 1) clustering within C2A1 and haplotype C5A4 (n = 1) clustering within C5A6 (Fig. 1).In both the C1A1/C2A1 and the C5A4/C5A6 haplotype clusters, the haplotypes were separated by only one base pair difference.Due to the proximity of these haplotypes, PERMANOVA analyses by haplotype were run with haplotype C1A1 merged with C2A1 and haplotype C5A4 merged with C5A6, as well as with haplotypes C1A1 and C5A4 excluded.In both analyses, there was a significant difference in microbiome composition among the haplotype groups (Merged: Pseudo-F, 3.7999, P = 0.0001; Excluded: Pseudo-F = 3.4659, P = 0.0001; Table S1).The C1A1/C2A1 cluster aligned with previous species predictions, as both haplotypes were assigned to Species Group 1 (Swierts et al. 2017).However, in the case of the C5A4/ C5A6 cluster, the haplotypes were assigned to Species Groups 3 and 4, respectively.Species Group 3 also contained haplotype C6A2, which in our case had a microbial community distinct from the other haplotypes.However, PER-MANOVA by Species Group was also significant (Pseudo-F = 3.0485, P = 0.001; Table S2), indicating that increased replication is required to resolve the placement of C5A4 into a Species Group with regard to its microbiome composition.Therefore, we excluded this sample from the subsequent analyses of microbiome differentiation found within the 4 Species Groups.
The most abundant microbial phyla across all the X. testudinaria samples were Chloroflexi, Crenarchaeota, Acidobacteriota, and Actinobacteriota (Fig. 2).Cyanobacteria were also abundant, although Species Group 3 had a lower relative abundance of Cyanobacteria compared to the other Species Groups.The most abundant class found in all Species Groups was Nitrososphaeria, within the phylum Crenarchaeota (Table S3).Crenarchaeota has previously been characterized as an important component in nitrogen cycling within barrel sponges (López-Legentil et al. 2010).Although previously, another nitrifying group of bacteria, Nitrospira, was identified as abundant within the Singapore giant barrel sponges (Swierts et al. 2018).Other abundant classes in the barrel sponges include Anaerolineae, Dehalococcoidia, TK30, and TK17 within the phylum Chloroflexi, Thermoanaerobaculia and Vicinamibacteria within phylum Acidobacteriota, Acidimicrobiia within phylum Actinobacteriota, Cyanobacteriia within phylum Cyanobacteria, BD2-11 terrestrial group within phylum Gemmatimonadota, and Gammaproteobacteria within phylum Proteobacteria.There were no differences in alpha diversity among Species Groups (P > 0.05; Table S4).We identified 52 core ASVs, when defined as presence in all X. testudinaria samples (Table S5), reduced from the 71 core ASVs found in 2018 (Swierts et al. 2018).Together the core ASVs represent 33.0% of the total sequence reads, with the most abundant phyla being Acidobacteriota (11.3%),Chloroflexi (9.28%), and Actinobacteriota (5.5%).
The GLM identified 85 ASVs, consisting of 12 Archaea and 73 Bacteria, as highly significant (P < 0.01) in differentiating among the four Species Groups (Table S6).The majority of significant ASVs were from the phylum Chloroflexi, of which there was a higher abundance in Species Groups 1 and 3 (Fig. 3).Species Group 1 also contained a higher abundance of significant ASVs from the following Phyla present in a relative abundance < 0.01 have been grouped together as rare taxa Fig. 3 The ASVs identified as highly significant (P < 0.01) in contributing to the differences in sponge Species Groups from the GLM, organized by phylum to demonstrate the relative abundance of the differentiating taxa for each species group (see Table S6 for detailed ASV taxonomy) 166 Page 6 of 8 phyla: Acidobacteriota, Bacteroidota, Gemmatimonadota, Nitrospirota, and Planctomycetota.Species Group 2 had a higher relative abundance of Archaea ASVs from the phyla Crenarchaeota and Thermoplasmatota.Species Group 3 had a higher relative abundance of ASVs from the phyla Actinobacteriota, PAUC34f, and Proteobacteria; while Species Group 4 had a higher relative abundance of ASVs from phylum Crenarchaeota and phylum Cyanobacteria.
Overall, these microbiome data support the placement of X. testudinaria into the Species Groups proposed by Swierts et al. (2017).Geographic differences in microbiomes exist across the sponges' broad range throughout the Indo-Pacific (Swierts et al. 2018), within a given region the Species Group differences in microbiome composition become evident.We recommend additional studies with increased sampling effort to include sufficient replication of all haplotypes within each geographic region to conclusively confirm the results reported here.Additionally, research coupling demographic surveys of X. testudinaria within Singapore with their genetic composition is required to determine if the barrel sponges are experiencing similar shifts in genetic structure as observed in the Caribbean (Deignan et al. 2018;Evans et al. 2021).These results highlight the usefulness of merging microbiome and molecular data in assigning species boundaries and hint at potential functional differences that allow the sponges to coexist sympatrically.

Table 1
Number of haplotypes found at each site, and the Species Group assignment of each haplotype according toSwierts et al.