Climate change threatens biodiversity worldwide, and assessing how those changes will impact communities will be critical for conservation. Dominant primary producers can alter local-scale environmental conditions, reducing temperature via shading and mitigating ocean acidification via photosynthesis, which could buffer communities from the impacts of climate change. We conducted two experiments on the coast of southeastern Alaska to assess the effects of a common seaweed species, Neorhodomela oregona, on temperature and pH in field tide pools and tide pool mesocosms. We found that N. oregona was numerically dominant in this system, covering >60% of habitable space in the pools and accounting for >40% of live cover. However, while N. oregona had a density-dependent effect on pH in isolated mesocosms, we did not find a consistent effect of N. oregona on either pH or water temperature in tide pools in the field. These results suggest that the amelioration of climate change impacts in immersed marine ecosystems by primary producers is not universal and likely depends on species’ functional attributes, including photosynthetic rate and physical structure, in addition to abundance or dominance.

Study Site: To evaluate the role of the abundant alga Neorhodomela oregona (Doty) Masuda in driving local climate conditions, we conducted “removal” and “mesocosm” experiments at John Brown’s Beach (57.05° N, 135.33° W) near Sitka, Alaska from 05 Jul 2019 to 27 Sep 2019. Southeast Alaska was an ideal location for this study as it has been subjected to relatively low levels of direct human disturbance yet is experiencing rapid environmental change (Stafford et al., 2000). Air temperature in Southeast Alaska has increased by ~0.11°C per decade since 1830 (Wendler et al., 2016; Jewett and Romanou, 2017), well above the global mean rate of 0.07°C per decade (since 1880; Blunden and Arndt, 2019). Sea surface ocean pH has declined by 0.03 units over a recent 15-year window (1991-2006) in the northeast Pacific waters off the Alaskan coast (Byrne et al., 2010).

Removal Experiment: We selected 10 tide pools, which ranged from 2.5 - 23.5 L in volume and 2.49 - 3.29 m in tide height (i.e., vertical position within the intertidal zone), which were separated by an average distance of 4 m, for the removal experiment. We began by assessing the physical characteristics of the experimental tide pools. We measured volume by pumping the water from a tide pool into a graduated bucket, and we assessed the basal surface area of the pool, as well as N. oregona abundance, by placing a flexible mesh quadrat with 10 cm × 10 cm squares on the bottom of each tide pool (Bracken and Nielsen, 2004; Sorte and Bracken, 2015; Silbiger and Sorte, 2018). Tide heights (in meters above mean lower low water) for each pool were measured using a sight level, a surveying rod, and tidal predictions from the National Oceanic and Atmospheric Administration (2019). We assigned pools to treatment and control groups (n = 5, removal or control) by randomizing assignments until various physical and biological metrics did not vary between treatment and control (based on a generalized linear model with a threshold of p > 0.2). Metrics included tide height, volume, basal surface area, percent cover of N. oregona, and species richness. We removed N. oregona from the treatment pools using scissors and by cutting as close to the substratum as possible without removing the holdfasts to avoid damaging surrounding organisms. We measured the wet biomass of N. oregona from each removal pool in the field before using that algal biomass in the mesocosm experiment (described below).

To assess the abundance of N. oregona and community composition in the tide pools, we conducted biodiversity surveys in the pools before and immediately after N. oregona removal (06 Jul – 19 Jul 2019), and then every two weeks until 27 Sep 2019 (for a total of seven surveys; Figures 1B & A1). During the surveys, we pumped water out of each tide pool, laid down a flexible mesh quadrat with 10 cm × 10 cm squares along the bottom, recorded the surface area covered by each sessile species (algae and invertebrates; 0.1 square or 10 cm² being the minimum measurement assigned for a species present in trace amounts) and counted all mobile invertebrates present (Bracken and Nielsen, 2004; Silbiger and Sorte, 2018). We identified organisms to the lowest possible taxonomic level: species when possible and genus when species were impossible to differentiate in the field (as with, e.g., Littorina plena and Littorina scutulata). In some cases, species were pooled and tallied together (e.g., “limpets” or “coralline algae”). 

To assess the impacts of N. oregona removal on tide pool pH, we conducted time-series samplings in the tide pools during the daytime and nighttime both before and after N. oregona removal (10 Jul – 16 Jul 2019; Figure A1). We measured temperature and salinity with a ProDSS Multiparameter Water Quality Meter (YSI, Yellow Springs, Ohio), and light intensity with a MQ-210 Underwater Quantum Meter (Apogee, Logan, Utah). Over the four sampling periods (day and night, both before and after removal), we took physical measurements at five time points over a ~2.5 h time series during low tide, sampling once every 30 min beginning immediately following isolation of the tide pools from the ocean. We also collected water samples at the first, third and fifth time points. The water samples were collected by hand-pumping 250 mL of water from the bottom of the tide pool into a vacuum flask, and then carefully siphoning the water into two 125 mL amber glass sample bottles to minimize gas exchange between the water sample and the atmosphere. All sample containers were rinsed three times with seawater prior to sample collection. We immediately preserved each water sample with 60 µL HgCl₂ and sealed them for later analyses to determine pH and total alkalinity.

We analyzed pH in the water samples from both experiments on a UV-1800 benchtop spectrophotometer (Shimadzu, Carlsbad, California), following best practices as described in Dickson et al. (2007). We divided each water sample into triplicate subsamples and analyzed them separately. We took an initial reading of each subsample at three wavelengths, added 50 µL of m-cresol dye, and mixed and re-analyzed the subsample at the same three wavelengths (Liu and Chan, 2011). We calculated the difference between the initial reading and the dye-added measurement, which we then used to calculate the pH value of each subsample. We took the mean of all subsamples with < 0.005 pH unit difference among them (subsamples outside that range were excluded) for each individual sample to produce a preliminary pH value. We then used CO2calc software (Robbins et al., 2010) to correct the preliminary pH value for total alkalinity (TA; analyzed as described below), salinity, temperature, and stoichiometric dissociation constants and calculate the final pH on the total scale (Mehrbach et al., 1973; Dickson and Millero, 1987; Kroeker et al., 2021).

We analyzed the TA of the water samples with open-cell titrations (as in Silbiger and Sorte, 2018) on a T50 titrator with LabX software (Mettler-Toledo AG, Schwerzenbach, Switzerland). We measured a certified reference material (Marine Physical Laboratory, Scripps Institution of Oceanography, La Jolla, California) at the beginning of each session as a standard (acceptable range: ±1% error), following an established protocol for open-cell TA analysis (SOP 3b) (Dickson et al., 2007; Silbiger and Sorte, 2018).

We conducted two additional samplings using a light-dark incubation method (Noël et al., 2010; Bracken et al., 2022) to assess how pH in the tide pools responded to differing light conditions. During these trials, we measured pH values across three time points using a HI9829 multiparameter meter with a 7609829 glass pH electrode (Hanna Instruments, Woonsocket, Rhode Island), which was calibrated using a Tris solution according to the best practices specified in SOP 6a by Dickson et al. (2007). We measured initial pH, remeasured following a ~30 min dark incubation period under an opaque, black plastic sheet, and collected a final measurement after a ~30 min light incubation period following the removal of the sheet.

Tide pool water temperatures were recorded every 5 min for the duration of the study by HOBO Pendant® Temperature/Light 64K Data Loggers (Onset Computer Corporation, Bourne, Massachusetts) anchored in the center of the pools. For comparison to our seawater temperature data, ambient air temperature data were sourced from the weather station at nearby Sitka Rocky Gutierrez Airport (Sitka, Alaska; < 1 km from the site) via CustomWeather, Inc. (2021).

 We conducted all statistical analyses in R (R-version 4.0.4; R Core Team, 2013) using a generalized linear mixed model (GLMM) repeated measures analyses and generalized linear models (GLM). We used a GLMM (‘lmer’ function; Bates et al., 2015) to evaluate the effect of the removal treatment on N. oregona abundance (cover) in the experimental tide pools and track recovery over time. N. oregona cover was modeled as a function of the fixed factors of treatment, time (bi-weekly surveys), and treatment × time, with the tide pool included as a random effect. We applied Kenward-Roger corrections to the GLMM to adjust the degrees of freedom to accurately reflect a repeated measures structure (Kenward and Roger, 1997; Kuznetsova et al., 2017) and conducted post hoc pairwise comparisons on N. oregona cover using Tukey’s HSD (‘emmeans’ function; Lenth, 2018).

To evaluate the effects of N. oregona on pH, we used the pH values at each of the three time points at which water was sampled to calculate the rate of pH change in tide pools (i.e., the slope of the relationship between pH and time), and we compared abundances of N. oregona to the calculated rate of pH change during the daytime and nighttime sampling periods. Similarly, we used the field pH measurements from the light-dark trials (which were subsequently converted from mV to pH units) to calculate the rate of pH change between the initial measurement and the measurement taken at the end of the dark incubation period to represent the rates of pH change during the night (Bracken et al., 2022), as well as the rate of pH change between the end of the dark incubation period and the final measurement (after a ~30 minute light incubation period) to correspond to the daytime water samplings. To assess the effects of N. oregona on water temperature in tide pools, we calculated the daily maximum water temperature for each tide pool over the full 11-week period following N. oregona removal.

We used GLMs (‘glm’ function in R) to assess the effects of N. oregona on pH. For intact pools prior to N. oregona removal, we evaluated the rate of pH change as a function of N. oregona area (in cm² of surface area per L of water volume), with the tide height of each pool, mean light in each pool (average of five time points; the light was not included in night analyses as it was uniformly measured as 0 at night), and mean water temperature in each pool during the sampling (across the five time points) included as covariates. Identical analyses were conducted on the pH data from the light-dark trials, with light intervals substituted for daytime samplings and dark intervals replacing nighttime samplings, except that individual temperature measurements were used rather than a mean value. This analysis of intact tide pools (before the removal) was also run with the assigned treatment group included as an additional factor, an analysis which confirmed that there was no initial difference in pH change between the treatment groups prior to removal (p > 0.4).

To test the effect of the N. oregona removal on pH, we evaluated the rate of pH change after removal as a function of treatment (removal vs. control), with tide height, mean water temperature, mean light, and pre-removal N. oregona area (in two-dimensional basal cover as measured in the biodiversity surveys) included as covariates, as well as an interaction between treatment and pre-removal N. oregona area. The interaction effect was included to assess whether the amount of N. oregona removed influenced the results, and we separately tested the effect of pre-removal N. oregona area in the removal and the control groups in the absence of other covariates to further investigate the role of initial N. oregona area as a potential driver of pH change. Finally, we conducted a combined analysis of the rates of pH change during day and night based on treatment, with pre-removal N. oregona area included as a covariate, as well as post-hoc tests comparing the treatment groups (‘emmeans’ function; Lenth, 2018). Assumptions of normality and homogeneity of variances were checked using the Shapiro-Wilk test and Levene’s tests, respectively.

We evaluated the role of the total producer and consumer assemblage in driving pH by comparing the pH change in each pool to total consumer abundance and producer dominance. Total consumer abundance was calculated using the surface area of all basal invertebrate species and converting counts of mobile invertebrates to surface area (Table A1). We did this conversion using photographic image analysis (with ImageJ; Abràmoff et al., 2004) of ~10 individuals per species of mobile invertebrate to find a mean surface area for an individual of each species and then multiplying that value by the number of individuals in each pool. For the few species we could not collect in the field, we substituted the measurements of species known to be of similar size (Table A1). We used 10 cm² as a minimum surface area for any mobile invertebrate species present, consistent with our methods used for the basal species in our community surveys. We then calculated consumer abundance as the total area per tide pool volume of non-photosynthetic species. Producer dominance, a metric used to represent the relative abundance of producers and consumers in an ecosystem, was calculated as the total abundance of all producer species (in two-dimensional basal cover from the biodiversity surveys) minus the total abundance of all consumers present. We modeled the rate of pH change as a function of total consumer abundance (cm² L^-1; ‘glm’ function) with tide height, mean water temperature, and mean light included as covariates, and ran similar analyses (with the same covariates included) on pH and producer dominance. Additionally, to account for the potential effects of the highly productive producer Ulva spp. (Sand-Jensen, 1988; Israel et al., 1995), we also ran the pre-removal and post-removal analyses of N. oregona abundance and pH with Ulva spp. abundance included as an additional covariate. The GLMs used in the removal experiment used a Gaussian distribution (identity link) except for the models of total consumer abundance and nighttime pH, which used a Gaussian distribution with an inverse link after the model failed to pass the Shapiro-Wilk test using an identity link.

To evaluate the effect of N. oregona removal on tide pool water temperature, we conducted a repeated measures analysis using a GLMM (‘lmer’ function, with Kenward-Roger corrections applied; Kenward and Roger, 1997; Bates et al., 2015; Kuznetsova et al., 2017) with data from the first month (prior to significant N. oregona recovery following the removal treatment; Figure 1) and, in a separate analysis, for the full 11-week duration of the study. The temperature was modeled as a function of the fixed factors of treatment, time (days), ambient air temperature, and an interaction between treatment and time, with the tide pool included as a random effect.

Mesocosm Experiment: We set up mesocosms on the beach adjacent to the experimental pools at our John Brown’s Beach study site on 13 Jul 2019. Mesocosms (12-L plastic tubs, n = 5 N. oregona addition, and n = 3 control) were arrayed in two parallel lines of four, randomly arranged with regards to treatment. We added N. oregona biomass from one of the n = 5 removal tide pools to each of the n = 5 addition treatment mesocosms. Each mesocosm also contained the quantity of seawater equal to the volume of the pool from which the N. oregona was removed (except that 10 L of seawater was added to the two mesocosms corresponding to the removal pools with >10 L volume). We added 10 L of seawater but no N. oregona biomass to the control mesocosms.

We conducted water sampling using a time series similar to the removal experiment (as described above), except that there was no “before” sample collection. We sampled the mesocosms after N. oregona addition during the daytime (4 h after algae were added to the mesocosms) and nighttime (10 h after addition; Figure A1). Prior to each time-series sampling, we simulated tidal inundation by flushing the mesocosms with seawater. We secured the algae in the mesocosms with wire mesh, poured the water out of the mesocosms, and used a graduated bucket to refill the mesocosms with the assigned volume of seawater. We took physical measurements at five time points over a ~2.5 h time series, sampling once every 30 min, and collected water samples on the first, third and fifth time points for later pH and TA analyses.

To test the effect of N. oregona on the rate of pH change in isolation, we applied GLMs (‘glm’ function) to the data from the mesocosms, for which we used two metrics of N. oregona abundance: source pool N. oregona surface area per mesocosm volume (cm² L^-1), which was the same metric we used for the algae in the field tide pools, and N. oregona biomass per mesocosm water volume (g L^-1), values that were only available for the mesocosms populated with the detached algae. We included mean water temperature as a covariate. We also used two GLMs (‘glm’ function) to analyze the combined day and night rates of pH change by treatment, with N. oregona biomass or source pool surface area of N. oregona included as a covariate, as well as post-hoc tests comparing the treatments in each model. Light measurements were not available for these analyses; however, the mesocosms were situated in an area of the beach with relatively homogenous light conditions (S. Mahanes, pers. obs.). Assumptions of normality and homogeneity of variances were checked using the Shapiro-Wilk test and Levene’s tests, respectively. All GLMs for the mesocosm pH analyses used a Gaussian distribution (identity link) except the analyses on the daytime sampling using biomass, which used a gamma distribution (inverse link) after the model failed to pass the Shapiro-Wilk test using a Gaussian distribution.