Disturbances can alter the structure and function of ecosystems. In stream ecosystems, changes in discharge and physicochemistry at short, intermediate, and long recurrence intervals can affect food webs and ecosystem processes. In this paper, we compare pH regimes in streams at La Selva Biological Station, Costa Rica, where episodic acidification frequency across the stream network varies widely due to buffering from inputs of bicarbonate-rich inter-basin groundwater. To examine the effects of acidification on ecosystem structure and function, we experimentally increased the buffering capacity of a headwater stream reach, and compared it to an unbuffered upstream reach. We compared these reaches to a naturally buffered and unbuffered reaches of a second headwater stream. We quantified ecosystem structural (macroinvertebrate assemblages on leaf litter and coarse woody debris), and functional responses (leaf litter and coarse woody debris decomposition rates, and growth rates of a focal insect taxon (Diptera: Chironomidae)). Non-metric multidimensional scaling and analysis of similarity revealed that macroinvertebrate assemblages were relatively homogenous across the four study reaches, although the naturally buffered reach was the most dissimilar. Ecosystem function, as measured by chironomid growth rates, was greater in the naturally buffered reach, while decomposition rates did not differ across the four reaches. Our results indicate that biological assemblages are adapted to pH regimes of frequently acidified stream reaches. Our experiment informs the effects on structure and function at short time scales in streams that experience moderate acidification, but larger magnitude acidification events in response to hydroclimatic change, as projected under climate change scenarios, may induce stronger responses in streams.

Study site

This study was conducted at La Selva Biological Station (LSBS, 10° 26’ N, 84° 01’ W), Costa Rica (Figure 1). La Selva (1536 ha area) consists of predominantly primary and secondary tropical wet forest, and receives, on average, 4000 mm of rainfall per year (Sanford et al. 1994). Streams at LSBS can be distinguished by two sources of groundwater inputs. Local groundwater inputs, which originate within the streams’ watershed, are present in all streams. Certain streams below ~50 meters above sea level receive inputs of inter-basin modified groundwater (IMG) (Pringle and Triska 1991, Pringle et al. 1993, Genereux et al. 2005). Streams with IMG inputs are characterized by higher solute concentrations, including soluble reactive phosphorus (SRP), cations, and carbonate species (Pringle et al. 1990, Pringle 1991, Oviedo‐Vargas et al. 2015). Streams with IMG exhibit faster microbial respiration (Rosemond et al. 2002, Ramírez et al. 2003), faster macroinvertebrate turnover (Ramírez and Pringle 2006), faster leaf litter decomposition rates (Ardón et al. 2006), and higher and more stable pH as a result of higher buffering capacity (Small et al. 2012).

Experimental design

We selected two headwater streams to test the effects of episodic acidification on both structure and function: Carapa and Arboleda Seep (hereafter ArbSeep) (Table 1). Carapa is a typical low-solute, headwater stream at LSBS, experiencing seasonal and episodic acidification (Small et al. 2012). ArbSeep is a headwater tributary to the Arboleda, and the downstream reach of ArbSeep receives IMG, resulting in different physicochemical characteristics between the upstream and downstream sections of ArbSeep. Due to logistical constraints, our replication was limited to two streams. We selected two 5-10 m reaches within each of the two streams, resulting in a total of four study reaches: two upstream low buffering capacity and frequently acidified reaches, one downstream high buffering capacity reach in ArbSeep (due to natural IMG inputs), and one downstream low-buffering capacity reach in Carapa that we chose for experimental carbonate amendments to simulate IMG inputs (Figure 1). The reaches in ArbSeep were separated by ~100 m, over which IMG inputs modify stream chemistry. In Carapa, study reaches were separated by ~15 m to aid in logistics in the field. The Carapa experiment occurred from 15 June 2018 – 20 July 2018, and measurements in ArbSeep occurred from 18 June 2018 – 23 July 2018. Rainfall data was collected from this period using a tipping bucket rain gauge located ~1 km from our streams and are available from the Organization for Tropical Studies (OTS, https://tropicalstudies.org/)

Buffering capacity addition

To mimic IMG inputs and prevent episodic acidification from occurring in the experimental downstream reach of Carapa, we continuously added buffering solution in the form of carbonate species in two phases for five weeks (15 June 2018 – 20 July 2018). For the first two weeks of the experiment, we dissolved 1,920 g of NaHCO₃ in a 220 L barrel of stream water and pumped this solution (~0.1 M HCO₃^-) into the stream above the downstream reach at a rate of 25 mL min^-1, refilling and amending the barrel with upstream water as needed. The NaHCO₃ pumping addition was not as effective as anticipated in preventing episodic acidification; therefore, for the remainder of the experiment, we added 500 g of solid CaCO₃ in coarse mesh bags directly in the stream immediately upstream of the downstream reach so it would dissolve in water. CaCO₃ bags were replaced every 2-3 days for the remaining three weeks of the experiment.

Stream physicochemistry

We deployed a YSI EXO1 multiprobe sonde (Xylem Inc., Yellow Springs, OH) in the two downstream reaches during the five-week experiment. Sondes measured temperature (°C) and pH every 15 minutes and were calibrated weekly. Due to sensor error, we only collected data in Carapa for the final three weeks of the experiment. As the sondes were deployed in the downstream reach during the experiment, we deployed the sondes in the upstream reaches for five weeks the following summer to collect the same data from the upstream reaches from August 13 - September 24 2019, during which there was 514.3 mm of rainfall. While these dates do not overlap with the dates of the 2018 experiment, the 2019 dates fall within the wet season at LSBS; in 2019, 811 mm of rainfall fell from June 15 - July 23 (Figure 2, a) Figure, while in 2018, over 1100 mm of rainfall fell during the measurement period (Figure 2, b). To examine the frequency of acidification, we used the Shapiro-Wilk test of normality to assess the distributions of pH data in each of the four study reaches, assuming that buffered reaches would exhibit normal distributions and upstream reaches would not. The STREAMS project at LSBS (a long term monitoring and biological research project) shows higher pH due to IMG inputs and skewed distributions of pH in streams with local groundwater only (Small et al. 2012, Ardón et al. 2013; Figure 1, b, c).

We measured discharge and dissolved nutrients weekly during the experiment. Discharge was measured in the downstream reaches using the cross-sectional area method with a velocity meter (Marsh-McBirney, Frederick, MA). Triplicate nutrient (NO₃^--N, NH₄⁺-N, and SRP) samples were collected in 60 mL bottles, filtered in the lab (0.45 µm), and frozen. Samples were transported to North Carolina State University frozen, where NH₄⁺ and SRP were analyzed on an AA3 Segmented Flow Analyzer (Seal Analytical, Mequon, WI) and NO₃^- was analyzed on a Metrohm 930 ion chromatograph (Metrohm, Ionenstrasse, Switzerland).

Individual organism response

We measured the growth of chironomid larvae (Diptera: Chironomidae) in the four reaches of the buffering experiment. Chironomids are the dominant detritivores in streams at La Selva and have comprised >50% of invertebrate biomass in previous leaf litter decomposition experiments (Rosemond et al. 2001, Ardón et al. 2006). Coarse mesh bags filled with air-dried Ficus insipida Willd. leaves were incubated in each of the four reaches for 10 days, allowing for colonization by chironomids, and returned to the lab. Chironomid larvae were picked from the leaf packs and measured to the nearest 0.2 mm. We selected six larvae between 2 and 5 mm for growth experiments. Groups of six larvae were placed into six replicate plastic tea strainers (mesh size = 224 µm) with six Ficus leaf discs (10 mm diameter) and deployed into the respective stream reach for 48 hr. After deployment, we counted surviving individual chironomids and re-measured the length of each individual. We estimated chironomid biomass using length-mass relationships (Small et al. 2011) as the mean of the initial six individuals deployed and the mean of the surviving individuals as the final biomass, and calculated percent change in biomass and survivorship for each tea strainer using the mean of the surviving individual chironomids.

Macroinvertebrate assemblage and ecosystem function response

In the four stream reaches, we deployed coarse mesh (3 mm) and fine mesh (300 µm) litter bags to measure leaf litter (LL) and coarse woody debris (CWD) decomposition. Fine mesh litter bags were used to exclude macroinvertebrates and estimate microbial decomposition rates, whereas coarse mesh bags allow macroinvertebrates access and contribute to decomposition (Woodward et al. 2012). Fine mesh bags were largely ineffective in excluding macroinvertebrates but were included in the analysis to determine what attributable effects remained. Litter bags of each mesh size and organic matter type were deployed for each collection date and type of organic matter, and triplicate litter bags were collected following one and five weeks of deployment in the stream (total n of decomposition bags = 96). Leaves from Luehea seemannii Triana & Planch. (Tiliaceae) (5 g) were collected as recently abscised litter and air dried for at least 3 days; CWD segments from Pentaclethra macroloba (Willd.) Kuntze (Fabaceae) (50 g) were collected from recently fallen trees and air dried for at least 3 days. In an attempt to standardize organic matter inputs across all decomposition bags, leaves of similar size and small sticks (circumference < 10 cm, length < 25 cm) were added to each bag. On the pre-assigned collection date, including a replicate set of sacrificed leaves and wood for each reach and mesh size on the deployment date, litter bags were removed from the stream using a dip net, transferred to a 4 L plastic bag, and transported to the lab. In the lab, leaf and wood matter were washed over a 250 µm sieve to remove all sediment and biota. Leaf litter and CWD material were placed into paper bags, dried at 40 ˚C for 48 h, and subsamples placed in a muffle furnace at 500 ˚C for 1 h to calculate dry mass (DM) and ash-free dry mass (AFDM), respectively. We calculated the decomposition rate constant for LL and CWD, k_LL or k_CWD, using a negative exponential model on AFDM data. Material collected in the sieve was preserved in 95% ethanol for macroinvertebrate identification. Macroinvertebrates were identified to the family level, where possible.

Statistics

Chironomid survival and growth rates were evaluated across the four study reaches. We used a nested ANOVA to analyze percent survival, percent biomass changes, and final chironomid biomass to determine differences between stream reaches. Percent biomass change was non-normally distributed and was log₁₀ transformed.

Decomposition rates for LL and CWD were analyzed using two-way nested ANOVA, which allows evaluation of differences between the two study streams, the two reaches in each stream, and the effect of bag mesh size. We used a post-hoc Tukey test to determine significant groups among the four reaches. Macroinvertebrate abundance from both fine and coarse mesh litter bags was measured after one and five weeks of deployment from the four reaches. Abundance data failed the normality assumptions for ANOVA; therefore we evaluated macroinvertebrate abundance across the four study reaches, mesh bag size, and organic matter type using the non-parametric Aligned Rank Transform test in the ARTool v0.10.8 (Kay and Wobbrock 2020) and post-hoc evaluation in the emmeans v1.5.3 R packages (Lenth 2020). To evaluate macroinvertebrate richness, we calculated the number of families in each litter bag and evaluated richness across the experimental factors using non-parametric Kruskal-Wallis rank sum test due to the non-normally distributed data.

To quantitatively assess the diversity of the macroinvertebrate assemblages over time, we used analysis of similarity (ANOSIM) using the Bray-Curtis dissimilarity index and 999 permutations for the treatments in the experiment. We used the ANOSIM R statistic to evaluate diversity across the experimental treatments (reaches, mesh size, organic matter), where R values closer to 0 indicate similarity versus values closer to 1 indicating dissimilarity. To document the diversity of the macroinvertebrate assemblage at the end of the experiment, we used non-metric multidimensional scaling (NMDS). In the NMDS, we used the Bray-Curtis dissimilarity index in three dimensions on log₁₀ transformed abundance data and visualized the ordination across our experimental treatments (study reaches, mesh bag size, and type of organic matter) to assess the influence of each treatment. Analysis of the macroinvertebrate assemblage was done in the vegan v2.5-7 R package (Oksanen et al. 2019).

We calculated the effect size statistic Hedges’ g to compare the upstream and downstream reach in each stream, isolating the effect within each stream of how the metrics in our experiment change moving from more frequently disturbed reaches into less frequently disturbed reaches. In ArbSeep, where no manipulations occurred, we used Hedges’ g to examine of the naturally buffered downstream reach relative to the unbuffered upstream reach. Similarly, in Carapa, we used Hedges’ g to examine the effect of the experimental buffering treatment in the downstream reach to the unbuffered upstream reach. We calculated g for metrics at the organismal (chironomid growth rate), community (macroinvertebrate abundance and richness), and ecosystem scales (k_LL and k_CWD). Hedge’s g is the preferred metric for small sample sizes (n < 20) and was calculated using the effsize v0.8.1 R package (Torchiano 2020).

The Rmd file contains all code to execute the analysis for this mansuscript, including pulling in data. Data are source from the same folder containing the Rmd file. The rendered file from the Rmd are also attached, both as .html and PDF files.