Dead coral skeleton (Siderastrea siderea) was collected from patch reefs in the upper Florida Keys using an underwater pneumatic drill with a 5 cm diameter core bit. Cores were returned to NOAA’s Atlantic Oceanographic and Meteorological Laboratory (AOML), cut into cylinders of consistent height (1 cm x 5 cm diameter), and treated with 12.5% sodium hypochlorite for 48 hrs to remove organic matter (Fang et al., 2013). Each was rinsed with deionized water, dried for 24 hrs at 60°C, dry weighed using a calibrated analytical balance (0.0001 g precision, Ohaus), and CT-scanned (SOMATOM, Siemens) to detect prior boring (Enochs et al., 2016a). Any skeletal cylinders with signs of prior bioeroding activity were not included in the study.

Fragments of reef rock (dead Siderastrea siderea) colonized by P. lampa (forma incrustans) and C. varians (forma incrustans) were collected from reefs in the Florida Keys (Cheeca Rocks, 24.8966 N, 80.6169 W) and State waters of Miami-Dade (Emerald Reef, 25.6742 N, 80.0987 W), respectively, using the same coring methodology. Samples were collected from seven P. lampa parent sponges (72 total samples) and six C. varians parent sponges (64 total samples). The number of samples collected from each parent sponge varied based on the size of the parent sponge, and ranged from eight to 16 samples. Care was taken to ensure each sponge fragment had at least one protruding oscula to conserve sponge filter-feeding activity. Fragments were transported to shore in thermally-insulated 75 L YETI^® coolers, cut to a consistent height (1.5 cm x 5 cm cylinder diameter; sponge tissue present at the base of the cylinder) using a diamond band-saw (Gryphon AquaSaw XL), and securely attached to the top of the skeletal cylinders using monofilament (0.8 mm diameter) and stainless-steel crimps ( Figure S1 ). The samples were then left in “recovery baskets” at the site of sponge collection for two-weeks to encourage healing of the scarred tissue and promote attachment to the skeletal cylinders ( Figure S1 ). Following the two-week recovery period, samples were returned to the University of Miami’s Cooperative Institute for Marine and Atmospheric Studies (CIMAS) and NOAA AOML’s Experimental Reef Lab (ERL).

Replicates were randomly allocated across eight independent aquarium systems (75 L glass tank circulating with 75 L sump tank) based on the maximum number of achievable replicates per parent sponge. Samples were acclimated to laboratory conditions for one week that were consistent with their in situ field conditions (28.1°C and 8.05 ± 0.02 pH; mean pH ± diurnal variability). Each tank received incoming seawater from Biscayne Bay (filtered to 1 μm) at a rate of 250 mL min^-1 to allow for complete system turnover every 15 hrs. Tank pH conditions were precisely controlled using two gas-specific (CO₂ and air) mass flow controllers (GFC series, Aalborg) in an automated system (Enochs et al., 2018). Each tank contained eight samples of C. varians, nine samples of P. lampa, and two control skeletal cylinders with no attached sponges that were used to assess passive dissolution. After acclimation, a one-week ramping period was used to bring the tanks to four randomly assigned pH treatments (two tanks per treatment): pre-industrial (8.15 ± 0.02; mean pH ± diurnal variability), contemporary (8.05 ± 0.02), and two future OA scenarios (moderate OA, 7.85 ± 0.02; extreme OA, 7.75 ± 0.02) (IPCC, 2019). Diurnal pH variability within each treatment mimicked that of natural reef ecosystems, with the pH maximum and minimum occurring at 6 pm and 6 am, respectively. Tanks were maintained at 27.5°C and supplied with 4.5 L hr^-1 of fresh seawater. Light was provided by 135 W LED arrays (Hydra 52 HD, Aqua Illumination), with peak photosynthetically active radiation (PAR) set to 250 µmol m^-2 s^-1. Sponges were fed a protein-rich macroalgal mix (N-rich High Pro, Reed Mariculture, Campbell, USA) every three hrs from 6 pm to 6 am using a custom doser. These daily increments supplied food to a final concentration of 25 μL L^-1 (Achlatis et al., 2017). Treatment conditions were maintained for 31 days.

Seawater pH (total scale) was measured continuously in each tank using low-drift Durafet pH sensors (Honeywell) that were calibrated twice weekly using spectrophotometric pH measurements (Agilent Cary, 8454 UV-Vis Spectrophotometer). Additionally, seawater samples (250 mL) were collected weekly from each tank for analysis of spectrophotometric pH, total alkalinity (TA; Apollo SciTech, AS-ALK2), and dissolved inorganic carbon (DIC; Apollo SciTech, AS-C3) in accordance with the manufacturers’ protocols (two seawater replicates; calibrated with certified reference materials, Scripps Institution of Oceanography). Parameters were used to calculate pCO₂ (μatm) and Ω_Ar with CO2SYS (Lewis and Wallace, 1998). Mean carbonate chemistry conditions across treatments and tanks are provided in Table 1 .

Buoyant weighing (Dodge et al., 1984) was used to quantify total bioerosion using a calibrated analytical balance (0.0001 g precision, Ohaus). Samples were placed on a suspended stainless-steel basket attached to the balance using tungsten wire (0.05 mm). Measurements were taken at the start of treatment conditions and 31 days after treatment conditions. Seawater density (Anton Paar, DMA 5000M) and temperature (Digi-Sense) were used to convert buoyant weight to skeletal mass. The diameter (cm) and height (cm) of sponge tissue was measured with calipers and used to calculate sponge tissue surface area (cm²) for standardization of the bioerosion rates.

Daytime and nighttime chemical dissolution was measured using the alkalinity anomaly technique (Zundelevich et al., 2007) 21 to 25 days into treatment conditions. Daily incubations were conducted using twenty custom-built light (i.e., clear acrylic) and dark (i.e., black acrylic) 750 mL chambers ( Figure S1 ). Each chamber contained a stir-bar magnetically attached to a 12 V motor that maintained continuous water motion. For each incubation period, four chambers (one sponge sample per chamber) were placed within a treatment tank (five tanks simultaneously used per incubation period), with the incubation order of tanks randomly selected. Sponge samples had a minimum of a one-day “rest period” between the day and night incubations. Sponge samples were monitored for signs of physiological stress (e.g., tissue paling, closing of oscula) after each incubation period. Dark incubations were run from 6:00 am to 9:00 am, with the treatment tanks blacked out the night before to ensure dark acclimation of the sponges prior to the incubations. Light incubations were run from 12 pm to 3 pm, with each chamber exposed to a consistent PAR (250 µmol m^-2 s^-1). Sponge samples were returned to their treatment tanks following each incubation, with the duration spent in the enclosed chamber never exceeding the 3 hr incubation period.

Chemical dissolution measurements were performed on a subset of each species (i.e., 48 per species, 12 per treatment, 6 per tank) from the 31-day experiment. Additionally, eight control skeletal cylinders were incubated to account for the impact of passive dissolution. Seawater samples (250 mL) were taken at the start and end of each incubation and analyzed for TA, pH, and DIC as above. Oxygen levels were monitored to confirm that incubations did not reach hyperoxic (200-250% dissolved oxygen; Cohen et al., 2017) or hypoxic (< 75% dissolved oxygen; Haas et al., 2014) conditions. Nutrients (ammonium, phosphate, nitrate, nitrite) were analyzed (SEAL Analytics) and applied as TA corrections to ensure that the change in TA used to estimate chemical dissolution came exclusively from calcium carbonate (CaCO₃) dissolution/precipitation and was independent of nutrient production/consumption (Webb et al., 2017). Biologically-mediated chemical dissolution (mg CaCO₃ cm^-2 d^-1) rates were calculated from the corrected change in TA using the formula described by Zundelevich et al. (2007) and were standardized to the tissue surface area of each sponge sample.

Daytime and nighttime mechanical bioerosion rates were evaluated through the collection of sponge sediment chips (Fang et al., 2013) 25 to 29 days into treatment conditions. The same subset of 48 sponges per species from the chemical dissolution experimentation were each placed within submerged, open-top containers (750 mL) for 10-hr day (7 am to 5 pm) and night (8 pm to 6 am) periods. Sponges had a one-day “rest period” between the day and night collections. Following the collection period, the seawater within each container was poured through a 100 µm sieve to remove coarser materials and vacuum filtered through a pre-combusted/pre-weighed GF/F glass microfiber filter paper (0.7 µm pore size). The filter paper was then combusted again at 550°C for 3 hrs, yielding only the weight of the collected sediment chips (Fang et al., 2013). Mechanical bioerosion rates (mg CaCO₃ cm^-2 d^-1) were standardized to sponge tissue surface area.

Prior to running the statistical tests, parent sponge identity and tank effects were evaluated for all bioerosion metrics. The assumptions of normality (Shapiro-Wilk test) and equal variance (Bartlett’s test) were met for both blocking factors, and statistical tests (ANOVA) determined that neither had a significant effect on any of the bioerosion metrics. This was further confirmed through additional testing with a reduced dataset of fully-crossed individuals across all tanks and treatments. Therefore, parent sponge identity and tank were removed as factors, with sponge sample used as a factor in the models. Separate one-way ANOVAs were used for the statistical analysis of bioerosion rates for each species with pH treatment used as a fixed effect. An individual ANOVA was conducted for all three bioerosion parameters: total bioerosion, chemical dissolution (day and night), and mechanical bioerosion (day and night) (α = 0.05). Post-hoc Tukey’s tests were used to evaluate treatment groups that were significantly different from one another (p < 0.05). Statistical analyses and development of plots were conducted with R (R Core Team, 2020), RStudio (RStudio Team, 2020), and ggplot extensions (Wickham, 2016).

Transcriptomic samples were taken from the 8.05 and 7.75 pH treatment groups at 12 pm on day five of the 31-day experiment. As to not affect downstream physiological measurements of the sampled sponges, small razor blade scrapes were taken, resulting in ~1-2 mm of sponge tissue. Prior to sampling, razor blades were cleaned using RNase AWAY™ (Thermo-Fisher Scientific). For transcriptomic profiling, only three of the seven P. lampa parent sponges and three of the six C. varians parent sponges were sampled. For each treatment (8.05 and 7.75 pH), this resulted in four replicates (two from each experimental tank) from each parent sponge. Tissue samples were rapidly transferred to 2-mL microcentrifuge tubes, submerged in 1 mL of TRIzol™ (Invitrogen), flash-frozen in liquid nitrogen, and stored at -80°C. Sampled sponges were monitored after the tissue collection, with full recovery (tissue healing) reported after two to three days.

Total RNA was extracted from each sponge biopsy following Mayfield et al. (2011), with an initial modification of a bead mill used in place of a mortar and pestle for sample homogenization. Frozen sponge tissues were thawed and stainless-steel beads (Navy Eppendorf RNA Lysis Kit) were added to each 2 mL sample tube. Homogenization was done using a FastPrep-24 Classic (MP Biomedicals) with four rounds of 30 s at 6 m/s. Between each FastPrep round, samples were placed on ice for 5 min to avoid heat buildup from homogenization. Following homogenization, 250 µL of homogenized sponge tissues and 750 µL of fresh TRIzol were combined in a new 1.5-mL microcentrifuge tube.

Tubes containing homogenized tissue and TRIzol were vortexed and incubated at room temperature (RT) for 10 min. Chloroform (200 µL) was added to each tube, and the samples were vigorously vortexed, incubated at RT for 3 min, and centrifuged at 12,000 xg and 4°C for 15 min. The RNA-containing aqueous phase (~550 µL) was transferred to new 1.5-mL tubes that were prefilled with 250 µL of “high salt” solution (0.8 M Na citrate and 1.2 M NaCl) and 250 µL of isopropanol. Tubes were vortexed and incubated on ice for 10 min. This was followed by centrifugation at 12,000 xg at 4°C (10 min each). The supernatant was removed, and the RNA pellet was washed with 1 mL of 75% ethanol. RNA pellets were dried on the lab bench for 60-100 min to evaporate any residual ethanol. Dried RNA pellets underwent a DNase I digestion with the following solution: 45 µL RNase-free water, 5 µL of 10x DNase I buffer (Machery-Nagel), and 0.5 µL DNase I (Machery-Nagel). RNA pellets were resuspended in 50.5 µL of DNase I solution and incubated at RT for 15 min.

The remaining RNA extraction steps were conducted using the RNA Clean & Concentrator™-5 kit (Zymo Research) with no deviations from the manufacturer protocol. RNA was eluted in 30 µL of nuclease-free water pre-warmed to 25°C. Quantification of RNA was assessed using a Nanodrop 8000 Spectrophotometer (Thermo-Fisher Scientific) and an initial quality check was made by electrophoresing 5 µL on 0.8% TBE-agarose gels with 1X SYBR^® Gold (Thermo-Fisher Scientific). Samples with > 25 ng/µL, and for which 28s and 18s bands were present on the gel, were shipped on dry ice to the Duke University Center for Genomic and Computational Biology. Additional quality control was undertaken by the sequencing core using a Qubit™ fluorometer (Thermo-Fisher Scientific) and an Agilent Bioanalyzer 2100 ( Table S1 ). Only samples with RNA integrity numbers (RIN) > 7.0 were prepared for RNA-Seq. cDNA libraries were prepped using Kapa Stranded mRNA-Seq library prep kit (Illumina) with 100 ng of total RNA. Samples were sequenced across one Illumina NovaSeq 6000 lane to produce 2 x 150 S-Prime base pair (bp) paired-end (PE) reads with requested 10 M PE reads per sample.

All scripts and data files used for the bioinformatic analysis are available on Github repositories for C. varians (Morris and Studivan, 2023) and P. lampa (Studivan, 2023). Demultiplexing was done by the sequencing facility, and read quality was assessed and visualized using FastQC (Andrews, 2010). Adaptor sequences and low-quality reads (Phred < 15) were removed using Cutadapt (Martin, 2011). Re-pairing and deduplicating of the trimmed reads was performed prior to the transcriptome assembly to remove any PCR duplicates. De novo transcriptome assembly was completed for both species using Trinity (Grabherr et al., 2011). The assembled transcriptomes were then filtered for reads ≥ 500 bp, mitochondrial RNA (mitoRNA) and ribosomal RNA (rRNA) contamination was removed, and GC-content was evaluated.

For C. varians, the sponge host and symbiont sequences were bioinformatically filtered and separated by comparison to the complete genomes of Amphimedon queenslandica (Demospongiae) and Symbiodinium microadriaticum (Symbiodiniaceae), as well as the NCBI’s non-redundant protein (nr) database. Annotation of the de novo assembled sponge transcriptomes was completed using UniProt (The UniProt Consortium, 2018), Gene Ontology (Ashburner et al., 2000), and KEGG (Kanehisa and Goto, 2000) annotation resources, and eggnog-mapper (Cantalapiedra et al., 2021) was used to generate annotations based on the GO/KOG databases.

Sponge reads were aligned and quantified to their respective de novo reference transcriptomes using bowtie2 (Langmead and Salzberg, 2012) and samtools before being read into R (R Core Team, 2020) and RStudio (RStudio Team, 2020) for downstream analyses.

An initial pre-filtering to remove low count genes (less than 1 count in greater than 8 samples) was used to reduce memory requirements and increase computational speed. Principal components analysis (PCA) was performed using a variance stabilized transformation (VST) (Love et al., 2014) of the low count filtered genes. The function PCAtools::eigencorplot (Blighe and Lun, 2020) was used to correlate metadata variables of interest with each PC axes using Pearson correlation coefficients (PCC) and statistical significance for the respective metadata variables. A likelihood ratio test (LRT) in DESeq2 (Love et al., 2014) was used to identify all differentially expressed genes between parent sponge for each species (p.adj < 0.05), with subsequent contrasts run between parent sponge (alpha < 0.05). To identify sets of unique and common genes between parent sponge contrasts within each species, Venn diagram analysis was performed using the R package venn (Dusa, 2022). To identify the effect of OA, the DESeq2 model ~ parent sponge + treatment was used, with genes deemed significantly differentially expressed if p.adj < 0.05. For all significance testing using DESeq2 (Love et al., 2014), the adjusted p-value was computed using the Benjamini-Hochberg methodology.

Visualization of the significant OA differentially expressed genes (DEGs) was done using the VST counts with parent sponge identity variance removed using limma (Ritchie et al., 2015). Variance removed VST counts were then used as input for visualization in ComplexHeatMap (Gu et al., 2016), with hierarchical clustering of genes (rows) and samples (columns) incorporated.

Gene ontology (GO) enrichment analyses were run using Cytoscape (Shannon et al., 2003) and BiNGO (Maere et al., 2005) using the gene to GO identifiers generated during the de novo transcriptome assembly annotations. To test for enrichment against the background gene set, the genes remaining after initial low count filtering were used. The hypergeometric test was utilized, and p-values were corrected with a Benjamini-Hochberg false discovery rate (FDR) with alpha set to 0.01.

Treatment had a significant effect on total bioerosion rates for zooxanthellate C. varians and azooxanthellate P. lampa (ANOVA: Table 2 , Figure 1 ). Post hoc analysis revealed that C. varians total bioerosion rates were significantly higher under 7.85 conditions than 8.15 (p = 0.012; Figure 1A ). While no other treatment groups were significantly different for C. varians, there was a trend towards increased bioerosion rates for sponges exposed to 7.85 pH conditions compared to 8.05 (p = 0.099). For P. lampa, total bioerosion rates measured in both OA treatments (7.75 pH and 7.85 pH) were significantly higher than the 8.15 treatment (p = 0.046 and p = 0.027, respectively; Figure 1B ). Relative to contemporary pH conditions, P. lampa bioerosion rates were elevated in the 7.75 and 7.85 OA treatments, although neither were statistically significant (p = 0.118 and p = 0.073, respectively).

Under 7.85 pH conditions, average rates were 1.03 ± 0.12 mg CaCO₃ cm^-2 d^-1 (mean ± SEM) for C. varians and 1.97 ± 0.31 mg CaCO₃ cm^-2 d^-1 for P. lampa. In comparison, a slight decline in bioerosion was measured for sponges in the most extreme 7.75 pH treatment group (0.85 ± 0.08 and 1.88 ± 0.24 mg CaCO₃ cm^-2 d^-1, respectively). The lowest average rates were found in the pre-industrial pH treatment (0.64 ± 0.05 and 1.05 ± 0.15 mg CaCO₃ cm^-2 d^-1, respectively). This was followed by the contemporary pH treatment group, where average bioerosion rates were 0.73 ± 0.07 mg CaCO₃ cm^-2 d^-1 for C. varians and 1.18 ± 0.13 mg CaCO₃ cm^-2 d^-1 for P. lampa.

The incubations of the control skeletal cylinders measured negligible passive dissolution (0.71 ± 0.31 μmol kg^-1; mean change in TA ± SEM) with no significant difference reported across treatments (F_3,4 = 1.47, p = 0.350). This TA signal represented only 7.9% and 6.1% of the average TA signal for C. varians and P. lampa, respectively. As such, chemical dissolution rates of the sponge samples were not corrected for passive dissolution.

Cliona varians daytime ( Figure 2A ) and nighttime ( Figure 2B ) chemical dissolution rates were not significantly different among treatments (ANOVA: Table 2 ). The highest average daytime rates were measured in the 7.85 and 8.05 treatments (0.038 ± 0.004 and 0.036 ± 0.005 mg CaCO₃ cm^-2 d^-1, respectively), whereas the lowest daytime rates were measured in the 7.75 and 8.15 treatments (0.024 ± 0.003 and 0.028 ± 0.004 mg CaCO₃ cm^-2 d^-1, respectively). In comparison, average nighttime chemical dissolution rates were higher under pre-industrial (0.031 ± 0.004 mg CaCO₃ cm^-2 d^-1) than contemporary pH conditions (0.022 ± 0.003 mg CaCO₃ cm^-2 d^-1), although the relationship was not significant. Sponges of the two OA treatments (7.85 and 7.75) were characterized by virtually identical nighttime chemical dissolution rates (0.028 ± 0.003 and 0.028 ± 0.002 mg CaCO₃ cm^-2 d^-1, respectively).

Pione lampa daytime chemical dissolution rates were significantly different across treatments ( Figure 2C ) whereas nighttime rates were not ( Figure 2D ) (ANOVA: Table 2 ). Post hoc analysis revealed that daytime chemical dissolution rates in the 7.85 treatment (0.054 ± 0.005 mg CaCO₃ cm^-2 d^-1) were significantly higher (p = 0.028) than those from the 8.15 treatment (0.029 ± 0.006 mg CaCO₃ cm^-2 d^-1). Average daytime chemical dissolution rates under 7.85 conditions were higher than that of the 8.05 (0.034 ± 0.006 mg CaCO₃ cm^-2 d^-1) and 7.75 treatments (0.034 ± 0.007 mg CaCO₃ cm^-2 d^-1), though neither were significant (p = 0.108 and p = 0.107, respectively). While nighttime chemical dissolution rates were not significantly different across treatments, sponges exposed to 7.85 pH conditions had higher average nighttime rates (0.057 ± 0.006 mg CaCO₃ cm^-2 d^-1) than those from the 8.15 (0.041 ± 0.006 mg CaCO₃ cm^-2 d^-1), 8.05 (0.039 ± 0.005 mg CaCO₃ cm^-2 d^-1), and 7.75 (0.042 ± 0.006 mg CaCO₃ cm^-2 d^-1) treatments.

Cliona varians daytime mechanical bioerosion rates ( Figure 3A ) were significantly different across treatments (ANOVA: Table 2 ). Daytime mechanical bioerosion rates under 8.05 pH conditions (0.052 ± 0.005 mg CaCO₃ cm^-2 d^-1) were significantly higher than those measured from the 8.15 treatment (0.022 ± 0.002 mg CaCO₃ cm^-2 d^-1) (p = 0.013). There was also a trend towards higher average daytime mechanical bioerosion rates in the 8.05 pH treatment than the 7.85 (0.041 ± 0.011 mg CaCO₃ cm^-2 d^-1) and 7.75 (0.037 ± 0.004 mg CaCO₃ cm^-2 d^-1) OA treatments, though neither difference was significant (p = 0.658 and p = 0.388, respectively). Average nighttime mechanical bioerosion rates ranged from 0.023 ± 0.002 mg CaCO₃ cm^-2 d^-1 (7.85 pH treatment) to 0.033 ± 0.004 mg CaCO₃ cm^-2 d^-1 (8.05 pH treatment) and were not significantly different across treatments (ANOVA: Table 2 , Figure 3B ).

Pione lampa daytime ( Figure 3C ) and nighttime ( Figure 3D ) mechanical bioerosion rates were not significantly different among treatments (ANOVA: Table 2 ). Average mechanical bioerosion rates were 0.091 ± 0.005 under daytime conditions and 0.089 ± 0.005 mg CaCO₃ cm^-2 d^-1 under nighttime conditions, and were similar across all treatments.

A total of 23 samples of C. varians, average sequencing depth of 10.2 x 10⁶ PE reads, and 22 samples of P. lampa, average sequencing depth of 11.6 x 10⁶ PE reads, were successfully sequenced. Trinity (Grabherr et al., 2011) de novo transcriptome assembly outputted a total of 346,222 and 220,699 sequences for C. varians and P. lampa, with an average contig length of 1,453 and 1,260 bp respectively ( Table 3 ). Quantification using bowtie2 (Langmead and Salzberg, 2012) resulted in counts for 45,009 and 88,405 genes for C. varians and P. lampa, respectively.

PCA identified parent sponge (PS) as the largest driver of gene expression, with a significant correlation of PS identity for C. varians on PC1 (r^2 = 0.66, explained variance = 69%) ( Figure 4A ) and P. lampa on PC2 (r^2 = 0.73, explained variance = 23%) ( Figure 4C ). For both species, visualization of PC1 (69% explained variance; 39% explained variance) and PC2 (11% explained variance; 23% explained variance) corroborated the correlation results, with no overlap between 95% confidence intervals for C. varians ( Figure 4B ) and P. lampa ( Figure 4D ). LRT analysis detected significantly differentially expressed genes (DEGs) when looking at PS identity for C. varians and P. lampa. Venn analysis identified 738 genes for C. varians and 963 genes for P. lampa that were present in all DESeq2 contrasts between PS ( Figure S2 ).

PCA revealed no strong correlation between pH treatment (8.05 pH and 7.75 pH) and gene expression for C. varians ( Figure 4A ) or P. lampa ( Figure 4C ). When accounting for PS identity, differential gene expression analysis identified sets of significantly DEGs between the 8.05 and 7.75 pH treatments for C. varians (177 DEGs; 9 upregulated and 168 downregulated) and P. lampa (105 DEGs; 17 upregulated and 88 downregulated) ( Table S2 ). Using generated annotation files, 95% (C. varians) and 77% (P. lampa) of significant DEGs received annotation (Supporting data, see methods section). Visualization of annotated genes identified clear treatment and expression clusters for C. varians ( Figure 5 ) and P. lampa ( Figure 6 ).

In C. varians, the most significantly upregulated DEG between the 8.05 and 7.75 pH treatments was carbonic anhydrase II [CA2, p.adj = 0.001, log2 fold change = 4.20]. Other significant upregulated DEGs included three genes involved in metal ion binding [cadmium/zinc-transporting ATPase HMA2, p.adj = 0.002, L2fc = 3.72; ulvan-active sulfatase, p.adj = 0.017, L2fc = 4.07; protein phosphatase 2C 8, p.adj = 0.011, L2fc = 3.00] and one gene linked to the production of heat shock proteins [heat shock cognate 71 kDA protein, p.adj = 0.001, L2fc = 2.75].

Significant downregulated DEGs between the 8.05 and 7.75 pH treatments, for C. varians, consisted of 25 genes related to transcription, including 23 genes involved in transcription initiation/regulation [e.g., nuclear transcription factor Y subunit alpha, p.adj = 0.031, L2fc = -4.88; DNA-directed RNA polymerase II subunit RPB1, p.adj = 0.015, L2fc = -5.43] and two genes specific to DNA binding [myb-like protein G, p.adj = 0.031, L2fc = -4.51; myocyte-specific enhancer factor 2C, p.adj = 0.016, L2fc = -3.74]. Eight downregulated DEGs were linked to the post-transcriptional modification of mRNA [e.g., E3 ubiquitin ligase SUD1, p.adj = 0.029, L2fc = -5.50]. A total of 26 downregulated DEGs were annotated to genes that are known to play a role in protein synthesis. This included 11 genes involved in translation initiation/regulation [e.g., translation initiation factor eIF-2b subunit gamma, p.adj = 0.045, L2fc = -4.99; eukaryotic translation initiation factor 3 subunit D, p.adj = 0.022, L2fc = -4.25], in addition to 15 genes associated with the production/maintenance of the 40S and 60S ribosomal subunits [e.g., 40S ribosomal protein S2, p.adj = 0.002, L2fc = -5.21; 60S ribosomal protein, p.adj = 0.016, L2fc = -6.50]. Other significant downregulated DEGs included 10 genes linked to protein degradation [e.g., polyubiquitin, p.adj = 0.026, L2fc = -1.69; ubiquitin-conjugating enzyme E2 J1, p.adj = 0.041, L2fc = -5.27], nine genes associated with cellular signal transduction [e.g., DEP domain-containing protein, p.adj = 0.013, L2fc = -5.48; CBL-interacting protein kinase 8, p.adj = 0.012, L2fc = -5.50], seven genes related to cell division [e.g., mitotic checkpoint protein bub3, p.adj = 0.032, L2fc = -3.82; cell division control protein 2 homolog A, p.adj = 0.043, L2fc = -3.23], four genes involved in sodium potassium ion transport [e.g., sodium/potassium-transporting ATPase subunit alpha-4, p.adj = 0.005, L2fc = -4.40], and three genes that mediate apoptosis [e.g., N-alpha-acetyltransferase 16, p.adj = 0.031, L2fc = -5.65; SNF-related serine/threonine-protein kinase, p.adj = 0.018, L2fc = -4.61].

For C. varians, GO enrichment of the significantly DEGs due to OA identified four significantly enriched GO terms all associated with the ribosome and ribosomal processes: Ribonucleoprotein complex (GO:0030529), Structural constituent of ribosome (GO:0003753), Cytosolic ribosome (GO:0022626), and Ribosome (GO:0005840). A list of genes associated with each term are provided in the Supporting Information (File S1).

In P. lampa, significant upregulated DEGs between the 8.05 and 7.75 pH treatments consisted of one gene involved with cellular signal transduction [TNF receptor-associated factor, p.adj = 0.043, L2fc = 0.36], one gene that plays a role in intracellular membrane trafficking and vesicular transport [Ras-related protein Rab-17, p.adj = 0.039, L2fc = 0.19], and two genes linked to protein degradation and protein turnover [BLOC-1-related complex subunit 8-like, p.adj = 0.043, L2fc = 0.18; ubiquitin, p.adj = 0.045, L2fc = 2.05].

The most abundant significantly downregulated DEGs for P. lampa included 12 genes related to the cell cycle, specifically 11 genes involved in cell division [e.g., G2 and S phase-expressed protein 1 isoform X4, p.adj = 0.001, L2fc = -0.48; wee1-like protein kinase 1-A, p.adj = 0.045, L2fc = -0.37], and one gene linked to cell differentiation [thioredoxin domain-containing protein 3 homolog, p.adj = 0.049, L2fc = -0.43]. Other downregulated DEGs consisted of nine genes that play a role in the production/regulation of microtubule motor proteins [e.g., tektin-1-like, p.adj = 0.041, L2fc = -0.37; coiled-coil domain-containing protein 170-like, p.adj = 0.048, L2fc = -0.43], five genes involved in the post-transcriptional modification of mRNA [e.g., histone pre-mRNA stem-loop binding, p.adj = 0.007, L2fc = -0.58; serine/arginine repetitive matrix protein 2-like isoform X1, p.adj < 0.001, L2fc = -0.35], five genes that mediate apoptosis [e.g., PTB domain-containing engulfment adapter protein 1, p.adj = 0.014, L2fc = -0.39; Villin-1, p.adj = 0.025, L2fc = -0.26], three genes associated with cellular signal transduction [e.g., Proline-rich transmembrane protein 1, p.adj = 0.017, L2fc = -0.16; LisH domain-containing protein ARMC9, p.adj = 0.017, L2fc = -0.44], and three genes involved in protein binding/transport [e.g., Kinesin-like protein KIF9, p.adj = 0.035, L2fc = -0.60]. There were an additional five genes involved in post-translational modification of proteins and protein turnover [e.g., GTP-binding protein 2, p.adj = 0.014, L2fc = -0.20; BRCA1 associated protein, p.adj = 0.046, L2fc = -0.26].

For P. lampa, GO enrichment analysis of the significantly DEGs due to OA identified only one enriched term: Axoneme (GO:0005930). Genes associated with this term are provided in the Supporting Information (File S1).

The present study identified a stimulating effect of OA on sponge bioerosion but indicate a potential threshold response that occurs regardless of the presence of sponge algal symbionts, as evident by an apparent nonlinear response seen for the zooxanthellate sponge, C. varians, and azooxanthellate sponge, P. lampa. These findings and measured rates are in support of prior sponge literature, with OA-accelerated bioerosion having been previously reported for zooxanthellate C. orientalis (Wisshak et al., 2012), C. celata (Wisshak et al., 2014), C. varians (Morris et al., 2022), C. caribbaea (Webb et al., 2017), and azooxanthellate C. delitrix (Morris et al., 2022) and P. lampa (Enochs et al., 2015). However, the majority of these studies detected a linear relationship between OA and sponge bioerosion (Wisshak et al., 2012, pH range 7.57 to 8.10; Wisshak et al., 2014, pH range 7.67 to 8.07), with only one (Enochs et al., 2015, pH range 7.74 to 7.98) describing a similar nonlinear relationship to our study (pH range 7.75 to 8.15). As we reported comparable results to Enochs et al. (2015) for P. lampa and C. varians, this could indicate that a nonlinear OA response is specific to these Caribbean sponge species, though further assessments are needed to determine whether this response is apparent for species in other regions, such as the Pacific.

Compared to total bioerosion rates, the nonlinear OA relationship was more pronounced for biologically-mediated chemical dissolution ( Figure 2 ), though this should be cautiously interpreted due to a lack of significance between the two OA treatments (7.75 and 7.85 pH). Contrary to prior studies that assessed the response for zooxanthellate Pacific species (Wisshak et al., 2012), OA-enhanced chemical dissolution was not evident for C. varians, with similar daytime rates measured in the moderate OA and contemporary pH treatments. Additionally, a trend towards decreasing chemical dissolution was observed under 7.75 pH conditions, as average rates from this treatment diminished to pre-industrial levels, which could represent bioerosion impairment under extreme OA conditions ( Figure 2A ). Nighttime chemical dissolution was considerably lower under contemporary pH conditions than the other three treatments ( Figure 2B ) and likely reflects a symbiont effect, as zooxanthellate species are known to have higher daytime chemical dissolution rates under present-day pH conditions (Fang et al., 2016; Webb et al., 2017; de Bakker et al., 2018). Webb et al. (2017) also found that relative to nighttime conditions, elevated daytime chemical dissolution rates are lost at high pCO₂ conditions, which they hypothesized as being due to an antagonistic relationship between enhanced photosynthetic efficiency of sponge photosymbionts under OA leading to a higher uptake of daytime CO₂. This in turn could counteract the benefits of acidified seawater conditions on daytime chemical bioerosion. While we did not measure photosynthetic efficiency in our study, this temporal shift across treatments was supported by our data, as similar daytime and nighttime chemical dissolution rates were observed for C. varians in the 7.75 pH treatment ( Figures 2A, B ).

The P. lampa chemical dissolution OA response differed slightly from that of C. varians, with markedly higher rates measured under moderate OA conditions (7.85 pH) than the three other treatments ( Figures 2C, D ), though this was only significant relative to the daytime pre-industrial rates. This implies that specific to P. lampa, the process of sponge chemical dissolution appears to be enhanced in the 7.85 pH treatment, but that a threshold is surpassed at lower pH conditions (7.75 pH) that curtails the positive influence of OA, in a similar manner as that of C. varians. This corroborates Enochs et al. (2015), where reduced chemical dissolution rates were measured for P. lampa in the lowest pH treatment (972.5 pCO₂) after ~22 hrs of treatment conditions, suggesting that the same nonlinear relationship to OA is apparent under short-term (i.e., < 1 day) and longer-term (i.e., > 21 days) OA exposure. Unlike C. varians, no temporal differences in daytime and nighttime chemical dissolution rates were detected for P. lampa, which supports prior characterizations of other Caribbean azooxanthellate sponge species (de Bakker et al., 2018).

Contrary to previous studies that described chemical dissolution as the primary sponge bioeroding mechanism sensitive to pH conditions (Wisshak et al., 2012; Wisshak et al., 2014), our data indicates that mechanical bioerosion, specifically under daytime conditions, may also be influenced by the pH of ambient seawater ( Figure 3 ). Contemporary daytime mechanical bioerosion rates were significantly higher for C. varians than the pre-industrial pH treatment ( Figure 3A ), an effect that was not detected for the azooxanthellate P. lampa ( Figure 3C ). This result was unexpected, considering that the process of chemical etching (i.e., chemical dissolution) of reef framework is thought to predate the mechanical removal of CaCO₃ (reviewed by Schönberg et al., 2017). Under this theory, accelerated chemical dissolution would likely increase the production of sponge carbonate chips, and therefore, lead to higher mechanical bioerosion rates. However, despite the lack of significance in chemical dissolution found across the pH treatments for C. varians, mechanical bioerosion rates were significantly higher under contemporary pH conditions. This suggests that seawater pumping rates responsible for the expulsion of sponge sediment chips could have been impacted by treatment conditions. While the underlying mechanisms behind this response need to be investigated further, this could be related to differences in the sponge energy budget across pH treatments (Fang et al., 2014). The energetically costly process of seawater pumping may be influenced by ambient seawater conditions, and therefore impact the rate at which carbonate chips are ejected from the sponge oscula. This may be further supported by the fact that only the zooxanthellate species showed this response in our study, as the transfer of autotrophic by-products are important components of the sponge energy budget (Weisz et al., 2010; Achlatis et al., 2019; Achlatis et al., 2021). A prior study found that prolonged OA exposure can significantly reduce the uptake of dissolved organic carbon by sponges (Achlatis et al., 2017), a result that could be indicative of reduced pumping activity under OA conditions (Achlatis et al., 2019). Alternatively, sponges could have retained the sediment chips under OA conditions (Emson, 1966) or reallocated resources to processes critical for survival, rather than bioerosion, as a result of physiological stress (Bates and Bell, 2017).

Short-term (five days) exposure to OA led to a change in gene expression of only a small proportion of genes, with 0.004% and 0.001% of the identified genes found to be significantly differentially expressed for C. varians and P. lampa, respectively. While the overall transcriptomic response was limited, the set of DEGs characterized for both species leads us to speculate that extreme OA conditions (7.75 pH) may be stressful for reef-excavating sponges. Prior studies have described the prominent downregulation of genes to be emblematic of a generalized stress response, as it is thought to be an adaptive mechanism for organisms to redirect energy resources towards maintaining critical physiological processes (Moya et al., 2012). This widespread downregulation has been reported for a non-bioeroding Caribbean Demospongiae (Histioneis tubifera) under temperature stress (Guzman and Conaco, 2016) and scleractinian corals (Acropora spp. and Pocillopora damicornis) under OA stress (Moya et al., 2012; Vidal-Dupiol et al., 2013). A similar stress response could be supported by our study, where 95% and 84% of all significant DEGs were downregulated for C. varians and P. lampa, respectively ( Table S2 ), though additional assessments with higher replication and deeper sequencing are needed.

The gene annotations for the significant downregulated DEGs further reinforce the proposed stress response to extreme OA. Significant DEGs implicated in sponge metabolism were prominently influenced by OA, with 19% (C. varians) and 34% (P. lampa) of all annotated downregulated DEGs in our study known to play a role in metabolic functioning. These included genes involved in carbohydrate, lipid, and amino acid metabolism, and suggests that sponges may suppress metabolism or modify metabolic activity in response to OA. The active suppression of metabolism has been described as a survival response to acute environmental stressors for coral adults (Moya et al., 2012), coral recruits (Albright et al., 2010), and sea-urchins (Todgham and Hofmann, 2009; O’Donnell et al., 2010), which could indicate that sponges may reallocate energy towards more important physiological needs in response to OA stress. Additionally, genes that participate in the negative regulation and mediation of apoptosis were significantly downregulated in the OA treatment for both species, which demonstrates that apoptotic processes increased under 7.75 pH conditions. DEGs specific to apoptosis have previously been observed for marine organisms in response to OA, thermal, and disease-induced stress (Libro et al., 2013; Guzman and Conaco, 2016; Dineshram et al., 2021).

While these downregulated DEGs may describe a generalized OA stress response for both species, there was a notable difference in the number of downregulated genes implicated in protein synthesis between species. For C. varians, protein synthesis was the biological function most impacted by OA, as these genes represented 29% of the significantly downregulated DEGs. This was corroborated by the GO enrichment analysis, as significantly enriched terms were linked to ribosomes, which are imperative for synthesizing proteins. In comparison, downregulated genes involved in protein synthesis were far less abundant for P. lampa, as they were only the sixth most impacted biological function, represented 7% of the significantly downregulated DEGs, and showed no significant GO enrichment for terms linked to protein synthesis. This could suggest that photosymbionts alter the molecular response of sponges to OA if they provide the sponge host with carbon-rich photosynthetic products (Achlatis et al., 2019), and/or impact the availability of ammonium (NH₄ ⁺) and bicarbonate (HCO₃ ^-), that may modify the sponge physiological needs resulting from OA stress. However, this could instead be broadly interpreted as response variability specific to each genus (i.e., Cliona vs Pione) and requires further evaluation.

Specific to the process of sponge chemical dissolution, the transcriptomic analysis identified a significant downregulation of DEGs involved in calcium:proton ion transport for C. varians (four genes) and P. lampa (four genes) under OA conditions, which is in support of a recently described OA response for C. varians (DeBiasse et al., 2022). Calcium ions (Ca²⁺) play an important role in chemical dissolution, as the transport of Ca²⁺ away from the dissolution site in exchange for protons (H⁺) in seawater is needed maintain an acidified microenvironment at the sponge/substrate interface to promote CaCO₃ dissolution (Sullivan et al., 1986; Webb et al., 2019). While Ca²⁺ is involved in other cellular roles in marine sponges, including tissue contraction and cell to cell communication (de Ceccatty, 1971; Lorenz et al., 1996), the downregulation of these DEGs could suggest that the intracellular exchange of Ca²⁺/H⁺ is less critical to the bioerosion process under OA, as the pH gradient between the sponge/substrate interface and ambient seawater is reduced. Alternatively, the downregulation of genes involved in Ca²⁺ transport could represent an overall decline in chemical dissolution in response to extreme OA conditions (at five days), as seen by the nonlinear bioerosion response reported in our study (at 21 to 25 days), though short-term chemical dissolution measurements are needed to test this hypothesis.

Of additional relevance to chemical dissolution, carbonic anhydrase 2 (CA2) was significantly upregulated in C. varians, with a fourfold increase in gene expression reported in the 7.75 pH treatment. This result is consistent with a previous study that described the upregulation of two carbonic anhydrase transcripts for C. varians (gamma stage) under short-term OA exposure (7.6 pH) (DeBiasse et al., 2022), indicating that this response is apparent across multiple C. varians growth forms (beta stage and gamma stage). Carbonic anhydrases (CAs) serve several important functional roles within sponges, one of which is catalyzing the interconversion of HCO₃ ^- and CO₂ involved in pH regulation at the site of carbonate dissolution (Webb et al., 2019). Carbonate ions (CO₃ ^-2) produced during chemical dissolution are converted to HCO₃ ^- and CO₂, with the ratio between the two depending on the rate at which protons (H⁺) are delivered and DIC is removed from the site of dissolution (Webb et al., 2019). CAs are thought to enhance the rate at which DIC is removed by promoting the protonation of HCO₃ ^- and therefore supporting the passive diffusion of CO₂ molecules away from the etching site (Webb et al., 2019). This in turn could accelerate chemical dissolution and may partially explain some of the underlying mechanisms behind OA-stimulated sponge bioerosion (Wisshak et al., 2012; Enochs et al., 2015).

However, despite this upregulation of CAs, chemical dissolution of C. varians was not significantly enhanced in the OA treatment. While a direct comparison between the transcriptomic and chemical dissolution measurements should be cautiously applied since they were derived from different time points (day five and days 21 to 25, respectively), this could indicate that CAs play less of a prominent role in chemical dissolution than previously reported. Alternatively, these findings could be the result of short-term upregulation of CAs that may have stimulated sponge bioerosion under OA at the time of transcriptomic tissue sampling. This effect may have diminished during the chemical dissolution measurements performed > 21 days into treatment conditions if the sponges were stressed from prolonged OA exposure. As these results represent a snapshot of the sponge transcriptomic response to OA at day five, further testing is needed to determine whether genes involved in chemical dissolution are differentially expressed under shorter or longer OA exposure to evaluate whether there is a significant correlation between CA gene expression and chemical dissolution. Additionally, since the upregulation of CAs was not evident for P. lampa, this instead could indicate that the role of CAs in sponge bioerosion may differ between zooxanthellate and azooxanthellate species. This could be supported by a prior study that measured 5-10x higher concentrations of carbonic anhydrase in zooxanthellate sponges compared to azooxanthellate species (Hill, 1996).

In our study, parent sponge was identified as a larger driver of sponge gene expression than pH treatment. This was not surprising, as several OA studies have previously described a similar genotypic effect. Glazier et al. (2020) found that genotype explained most of the gene expression variance for the scleractinian coral, Leptoseris pertusa, when exposed to 6 months of OA, a result that was supported by physiological data from the same experiment (Kurman et al., 2017). This type of molecular response has also been described for parental genotypes of spiny damselfish, Acanthochromis polyacanthus (Schunter et al., 2016). While our study is the first to characterize DEGs parent sponge variability for both zooxanthellate and azooxanthellate bioeroding sponges under OA, and therefore this variability needs to be investigated further, this could imply that some sponges are more resilient to OA, in a similar manner to that of thermally resilient coral genotypes (Barshis et al., 2013). Contrary to the studies by Glazier et al. (2020) and Kurman et al. (2017), however, we did not observe a corresponding significant effect of parent sponge identity on the physiology data, as sponge bioerosion rates were not significantly different between parent sponges. This suggests that OA-induced parent sponge variability in gene expression may not have a considerable influence on sponge bioerosion. However, this could also be related to the duration of our experiment, where 31 days of treatment conditions may not have been long enough to fully capture the physiological impact of this parent sponge variability.

Contrary to prior studies, these data provide a less-than-dire impact of OA on sponge bioerosion under projected end-of-the-century pH scenarios. While the buoyant weight-derived total bioerosion rates ( Figure 1 ) support previously reported OA-stimulated sponge bioerosion (Wisshak et al., 2012; Enochs et al., 2015; Morris et al., 2022), the chemical dissolution rates ( Figure 2 ) and transcriptomic profiling ( Figures 5 , 6 ) indicate that this enhancement could be unsustainable and may plateau or diminish in response to prolonged OA. The molecular data identified an onset of stress related DEGs for both species after only five-days of extreme OA exposure, suggesting that the sponges may have already experienced the effects of physiological impairment early in the experiment. This OA stress response was further reinforced by the chemical dissolution data, where a trend towards depressed bioerosion was apparent in the 7.75 pH treatment, though this should be cautiously interpreted due to a lack of significance compared to the 7.85 pH treatment. While reduced bioerosion in the extreme OA treatment relative to the moderate OA treatment was still present in the total bioerosion rates ( Figure 1 ), the drop-off was less notable than that of chemical dissolution, which is likely a function of differences in response data that these two methodologies are designed to address.

The buoyant weight methodology could be less effective at capturing a sponge stress response, since the bioerosion rates measured using this methodology were calculated exclusively from change in mass after 31 days of treatment conditions. Therefore, OA-depressed sponge bioerosion captured using this technique could have been muted if the sponges were able to maintain OA-stimulated bioerosion in the short-term, before eventually succumbing to OA-related stress further into the experiment. By contrast, the 3 hr incubations used to measure chemical dissolution were designed to capture instantaneous bioerosion rates and therefore may better reflect the physiological state of the sponges after 21 days of OA exposure (though food depletion, waste accumulation, and acidosis effects could have influenced the measured rates, see de Bakker et al., 2018). This aspect of the chemical dissolution methodology may also account for discrepancies between our findings and that of other sponge bioerosion studies, which used rapid techniques and relatively short exposure to assess responses to OA (4 days see Wisshak et al., 2014; 22 hrs see Enochs et al., 2015), suggesting that sponge bioerosion may not have been impaired by OA stress at that point in their experiments.

Based on these results, predictions calling for sponge bioerosion to double by the end of the century may be unrealistic if the sponges become physiologically impaired by extreme OA conditions (Enochs et al., 2015; Bates and Bell, 2017). This prognosis is further supported by previous studies that found a similar physiologically damaging effect (e.g., tissue necrosis, bleaching, mortality) from predicted “business-as-usual” seawater temperature conditions (Fang et al., 2014; Fang et al., 2018; Ramsby et al., 2018). However, while this may ameliorate reef erosion in the long-term, our data agrees with prior predictions for accelerated sponge bioerosion under moderate OA conditions (Enochs et al., 2015). When considering that OA is simultaneously expected to reduce coral calcification (Langdon and Atkinson, 2005; Chan and Connolly, 2012) and enhance the bioeroding activity of other reef organisms (Tribollet et al., 2009; DeCarlo et al., 2015; Enochs et al., 2016a), this accelerated sponge bioerosion under moderate OA conditions could still have considerable implications for the persistence of coral reef ecosystems in the upcoming decades, even if OA-enhanced bioerosion plateaus or diminishes later in the century.