Marine deforestation leads to widespread loss of ecosystem function

Trophic interactions can result in changes to the abundance and distribution of habitat-forming species that dramatically reduce ecosystem functioning. In the coastal zone of the Aleutian Archipelago, overgrazing by herbivorous sea urchins that began in the 1990s resulted in widespread deforestation of the region’s kelp forests, which led to lower macroalgal abundances and higher benthic irradiances. We examined how this deforestation impacted ecosystem function by comparing patterns of net ecosystem production (NEP), gross primary production (GPP), ecosystem respiration (Re), and the range between GPP and Re in remnant kelp forests, urchin barrens, and habitats that were in transition between the two habitat types at nine islands that spanned more than 1000 kilometers of the archipelago. Our results show that deforestation, on average, resulted in a 24% reduction in GPP, a 26% reduction in Re, and a 24% reduction in the range between GPP and Re. Further, the transition habitats were intermediate to the kelp forests and urchin barrens for these metrics. These opposing metabolic processes remained in balance; however, which resulted in little-to-no changes to NEP. These effects of deforestation on ecosystem productivity, however, were highly variable between years and among the study islands. In light of the worldwide declines in kelp forests observed in recent decades, our findings suggest that marine deforestation profoundly affects how coastal ecosystems function.


Introduction
Consumers fundamentally affect ecosystems through trophic interactions [1]. These interactions are especially important if they result in changes to the abundance or distribution of ecosystem engineers, such as forest-forming trees, which can lead to changes in microclimates, biodiversity, primary production, nutrient cycling, and energy flow [2]. For example, the reintroduction of gray wolves (Canis lupus) into Yellowstone National Park, USA in the 1990s resulted in increased predation on elk (Cervus elaphus) and subsequently reduced herbivory on canopy-forming trees such as aspens (Populus tremuloides), willows (Salix spp.), and cottonwoods (Populus spp.) [3]. This ultimately led to changes in the morphology and hydrology compensatory production by any remaining fleshy macroalgae, encrusting coralline algae, and microalgae [28][29][30], which can result in greater NEP. Thus, understanding how GPP, Re, and NEP change with kelp forest change can be instrumental in discerning the broader impacts of deforestation on ecosystem productivity. This may be especially relevant for the Aleutian Archipelago where widespread kelp deforestation has resulted in significant reductions in fishes, invertebrates and fleshy macroalgae, increases in the exposure of encrusting coralline algae [12,31], and elevated benthic irradiances [14].

Results
We used benthic chambers to study patterns of GPP, Re, and NEP within remnant kelp forests, urchin barrens, and habitats that were in transition to becoming urchin barrens (i.e., they had lost all benthic fleshy macroalgae but still had abundant stands of the canopy-forming Eualaria fistulosa; Fig 1) at nine islands spanning more than 1000 kilometers of the Aleutian Archipelago (Fig 2, Table 1). Kelp forests and urchin barrens occur as alternate stable states of one another, often with sharply delineated boundaries between them, and exhibit little-to-no overlap in community assemblages [15,33] (Fig 1). Indeed, the benthic communities within our chambers reflected these assemblages, with the chambers deployed in the kelp forests having more than a 10-fold greater biomass of fleshy macroalgae, which were predominantly stipitate kelps, than those deployed in the urchin barrens, and the chambers deployed in the urchin barrens having a nearly 3-fold greater biomass of urchins than those deployed in the kelp forests (Figs 3 and 4). The chambers deployed within transition habitats contained high abundances of urchins and little-to-no fleshy macroalgae, except for the canopy-forming E. fistulosa. The chambers within all three habitats had high bottom covers of encrusting coralline algae below the fleshy macroalgae, which became exposed following deforestation. Benthic irradiances, measured as photosynthetically active radiation (PAR), varied among the three habitat types (ANOVA: F 2,14 = 4.826, p = 0.025), but this was variable among the nine islands and two study years (Habitat � Island(Year) interaction: F 14,33 , = 4.426, p <0.001; Table 2). Generally, PAR was greatest in the urchin barrens, lowest in the kelp forests, and intermediate in the transition habitats ( Fig 5). We examined how GPP, Re and NEP, and the balance between GPP and Re differed among the habitat types by measuring changes in seawater oxygen concentrations within replicate (n = 3) chambers (collapsible benthic incubation tents; hereafter cBITs) that were placed on the benthos over representative assemblages within each habitat type at each island. We predicted that NEP at the benthos would be reduced in the urchin barrens due to the loss of photosynthetic macroalgae. Instead, we found that NEP did not differ between the habitat types ANOVA: F 2,14 = 0.530, p = 0.600), nor did it differ from zero (i.e., GPP = Re) in any of the Photographs of each habitat type showing (A) high abundance of benthic macroalgae and canopy-forming kelps in the kelp forests, (B) lack of benthic macroalgae but remaining canopy-forming kelps and high abundances of sea urchins in the transition habitats, and (C) lack of benthic macroalgae and canopy-forming kelps, but high abundances of sea urchins in the urchin barrens.  Table 2). Specifically, NEP was greater in the urchin barrens than in the kelp forests on five of the islands and lower in the urchin barrens on three of the islands, with the average difference being 25.58 ± 373.26 mg O 2 m -2 day -1 lower in the urchin barrens (i.e., the deforested habitats) (Tables 3 and 4). The change on one island (Attu) was not determined due to lost replication (Tables 1 and 4). However, when averaged across all nine islands, NEP was generally lowest (-239.73 ± 425.16 mg O 2 m -2 day -1 , mean ± SE) in the kelp forests, highest (-59.60 ± 145.32 mg O 2 m -2 day -1 ) in the transition habitats, and intermediate (-120.08 ± 338.07 mg O 2 m -2 day -1 ) in the urchin barrens (Table 3, Fig 4). Benthic GPP also did not vary among the habitat types (ANOVA: F 2, 14 = 0.234, p = 0.794), but when averaged across islands, GPP was highest in the kelp forests (1,806.14 ± 521.75 mg O 2 m -2 day -1 ; mean ± SE), lowest in the urchin barrens (1,367.77 ± 483.99 mg O 2 m -2 day -1 ), and intermediate in the transition habitats (1,494.22 ± 452.41 mg O 2 m -2 day -1 ) (Fig 4; Table 3). Like NEP, the effects of habitat type varied among the nine islands visited in the two study years (Habitat � Island(Year) interactions: F 14,41 = 2.166, p = 0.028; Table 2). Specifically, GPP was lower in the urchin barrens than in the kelp forests on all but two of the islands, by an average of 461.60 ± 578.69 mg O 2 m -2 day -1 (mean ± SE) ( Table 4). Re also did not vary among the habitat types (F 2,14 = 0.390, p = 0.684), but when averaged across all nine islands, Re was again highest in the kelp forests ( Table 2. Results of separate three-way nested analyses of variance testing for differences in A) net ecosystem production (NEP), B) gross primary production (GPP), C) ecosystem respiration (Re), D) the range between GPP and Re, and E) irradiance (PAR) among the two sample years, nine islands, and three habitat types (kelp forests, transition habitats, and urchin barrens). For each analysis, year and habitat type were fixed factors, and island nested within year was a random factor. The model r 2 is given for each analysis.  Table 3). As with NEP and GPP, the effects of habitat type varied among the nine islands visited in the two study years (Habitat � Island(Year) interactions: F 2,14 = 2.744, p = 0.006; Table 2). Specifically, Re was lower in the urchin barrens than in the kelp forests on four of the islands and greater in the urchin barrens on four of the islands, with the average difference being 472.09 ± 734.70 mg O 2 m -2 day -1 lower in the urchin barrens (Table 4). Lastly, the range between GPP and Re, which we believe to be a better measure of ecosystem function regarding productivity than NEP alone, did not differ among the habitat types (ANOVA: F 2,14 = 0.318, p = 0.733), but was again greatest in the kelp forests ( Table 2). Specifically, the range between GPP and Re was lower in the urchin barrens than in the kelp forests on five of the islands and greater in the urchin barrens on two of the islands, with the average difference being 933.69 ± 1,262.65 mg O 2 m -2 day -1 lower in the urchin barrens (Table 4).

A) NEP
Although the effects of deforestation on all three metrics of productivity varied among the islands visited in the two study years, some general patterns were evident. When considered across all nine islands, GPP, Re and the range between GPP and Re were each greatest in the kelp forests, intermediate in the transition habitats, and lowest in the urchin barrens. Specifically, GPP was 24% higher, on average, in the kelp forests than in the urchin barrens, and 17% higher, on average, in the kelp forests than in the transition habitats, but it differed by only 7% between the transition habitats and urchin barrens ( Table 3, Fig 4). Benthic Re was 26% higher, on average, in the kelp forests than in the urchin barrens, and 22% higher in the kelp forest than the transition habitats, but it differed by less than 1% between the transition habitats and  the urchin barrens. The range was between GPP and Re was 24% greater, on average, in the kelp forests than in the urchin barrens, and 19% greater in the kelp forests than in the transition habitats, but it varied by less than 6% between the transition habitats and the urchin barrens. In contrast, GPP was, on average, greatest in the urchin barrens, intermediate in the transition habitats, and lowest in the kelp forests in contrast, but it did not differ from zero (i.e. GPP = Re) in any of the habitats. These patterns, however, were highly variable among the different islands visited in the two study years for each of the production metrics. Altogether, this indicated deforestation resulted in widespread but geographically variable losses to primary production and respiration by the ecosystem.
As with previous studies in aquatic ecosystems, we found that GPP and Re are generally in balance, resulting in exhibit GPP / Re ratios near 1.0, and NEP values near zero [21,22]. When examined within each cBIT separately, GPP and Re were consistently similar in magnitude with no differences in GPP / Re ratios among habitat types (ANCOVA: F 2,62 = 0.16, p = 0.852) ( Table 5, Fig 6). Further, the distribution of these ratios was symmetrical around 1.0 in each habitat (Fig 7). Interestingly, the highest individual values of NEP were not observed in the kelp forests but rather in the urchin barrens, which we believe was due to higher irradiances in the urchin barrens than the other two habitats (Fig 5) combined with compensatory production by the encrusting coralline algae and benthic diatoms [30]. However, those few observations aside, it is clear that all three benthic habitats remain in balance following deforestation, with GPP � Re, GPP / Re ratios � 1, and median NEP values � 0. Thus, although NEP may help differentiate between productive and unproductive ecosystems [22], it poorly describes changes in primary productivity following large-scale habitat change in the Aleutian Archipelago. Instead, it is clear that deforestation results in significant changes to the region's benthic communities, and these led to geographically variable reductions in GPP, Re and the range between them, which better reflect a reduction in ecosystem functioning. Further, it appears that even partial deforestation, where the benthic macroalgae and invertebrates have been lost but the canopy-forming kelps remain, results in lower GPP and Re at the benthos that is similar to trends found in urchin barrens.

Discussion
Trophic interactions can lead to changes in the distribution and abundance of habitat-forming species, which can have profound impacts on ecosystem function [2,9]. Deforestation, in particular, can result in changes to biodiversity and energy flow [2], altered regional and global climates [34], and even lead to species extinctions [35]. Coastal kelps are a pertinent example of such ecosystem engineers in nearshore habitats that have suffered declines in some locations over the past few decades due to both biological and physical stressors [10,11,[36][37][38][39]. Consequently, while kelp populations remain stable in many of the world's ecoregions [10,40,41], or may even be expanding in some high latitude regions in response to ocean warming [39,42], our study is relevant to many areas of the world where kelp forests have exhibited local to broad scale declines [10,[43][44][45][46][47]. Indeed, recent estimates suggest that global declines in kelp abundances may be as high as 2% per year [11], which can negatively affect numerous other species that depend on them for food and habitat. Certainly, the kelp forests of the Aleutian Archipelago are in critical condition in the face of widespread overgrazing by urchins, and this  (Table 1). Each data point represents measurements from a single cBIT. Note the urchin barrens have the highest ratios observed, and the kelp forests have the largest number of low values. The vertical dashed line represents the 1:1 ratio. https://doi.org/10.1371/journal.pone.0226173.g007

Table 3. Community productivity values (measured as mg O 2 m -2 day -1 ) for A) Net ecosystem production (NEP), B) gross primary production (GPP), C) ecosystem respiration (Re), and D) the Range between GPP and Re (Range) estimated for each habitat on each island. Data reflect the means (SD) of the replicate chambers in each habitat (kelp Forests, Urchin Barrens and Transition Habitats). Positive values for NEP reflect net oxygen production and negative values reflect oxygen consumption (net respiration). Negative values for
Re reflect oxygen consumption (i.e. more negative values reflect greater respiration by the ecosystem). "NA" denotes not available due to lack of replication (i.e. data are based on only one chamber at that island; see Table 1).

Island Kelp Barren Transition
Adak -851. 56  has had profound effects on the region's benthic communities and on patterns of gross primary production and ecosystem respiration. Whether these forests will recover and return to prior ecosystem functioning regarding these metrics is unknown, but observations of kelp forests from other areas of the world suggest it is possible. For example, Laminaria longicruris forests recovered from overgrazing following localized disease outbreaks that decimated sea urchin populations in Nova Scotia [48], while L. hyperborea forests recovered in mid-Norway due to low sea urchin recruitment [49]. Ecklonia maxima expanded its range eastward in South Africa, coincident with cooling of the local ocean waters [50]. Likewise, Macrocystis pyrifera recovered along a~100 km stretch of the Pacific coast of Baja California, Mexico following nearly two decades of absence after the strong 1997-98 El Niño Southern Oscillation [51].
Recovery of the Eualaria fistulosa forests throughout the Aleutian Archipelago, however, would likely require widespread mortality in the urchin populations, which today seems unlikely. One potential contributing factor for this may lie in the low abundance of other urchin predators, such as the urchin eating sea star, Pycnopodia helianthoides [38,52,53], which historically has not been found in high abundances in the central or western Aleutians. Therefore, until predation on the urchins recovers or the urchin populations suffer widespread disease that reduces their numbers, benthic algal abundances, GPP and Re will likely remain generally lower in areas of kelp forest loss because the high abundance of urchins limits regrowth of macroalgae and maintains the urchin barrens [15]. Thus, we present a benchmark against which we can evaluate this recovery if it occurs, and understand the effects of further deforestation in this ecosystem. Although we have learned much about the effects of the otter-urchin-kelp trophic cascade in the Aleutian Archipelago, this study offers new insights into the consequences of such widespread deforestation on the region's benthic primary productivity. Certainly, benthic GPP, Re and the range between them are generally greatest in the kelp forests where macroalgae, fish, invertebrates, and microbial communities are all most abundant [15,[23][24][25][26]33], while they are lowest in the urchin barrens. Deforestation thus resulted in overall reductions in each of these metrics, identifying a general loss of ecosystem function. This, however, was geographically variable, with some islands showing elevated primary production following deforestation, which we believe is due to higher irradiances combined with compensatory production by microalgae (e.g. diatoms) and the coralline algal crusts. Indeed, we observed some of the highest production values in a few of the barrens cBITs where diatom mats formed within the chambers during the deployments. These cBITs also tended to have low numbers of urchins within them, and the chambers therefore appeared to exclude urchins from grazing the microalgae. In contrast, benthic primary productivity and respiration by the ecosystem are all similar in the urchin barrens and transition habitats, which have similarly high abundances of urchins and low biomasses of macroalgae [15,33], suggesting that the transition habitats have already suffered reduced ecosystem functioning. This, of course, reflects patterns at the benthos and

Table 4. The effects of habitat change on patterns of productivity for A) Net ecosystem production (NEP), B) gross primary production (GPP), C) ecosystem respiration (Re), and D the Range between GPP
and Re (i.e. Range) estimated for reach island. "Change" reflects absolute differences in each metric (measured as mg O 2 m -2 day -1 ) as the habitat transitions from Kelp forests to Transition Habitats, Transition Habitats to Urchin Barrens, Kelp Forests to Urchin Barrens (i.e. the total change due to deforestation). Positive values denote greater values for that metric and negative values denote lower values for that metric. NA denotes comparison "not available" due to loss of replicates in one habitat that precluded reliable estimates of the change (see Table 1). At the bottom of each not in the mid-water or at the surface where the canopy-forming Eualaria fistulosa remains abundant in the transition habitats. It is likely that these canopy-forming macroalgae would enhance GPP and perhaps result in positive values of NEP in the mid-water and at the surface in both the kelp forests and transition habitats. However, at the benthos, GPP and Re remain in balance following deforestation, leading to similar, near-zero NEP in all three habitats. We believe this reflects balance between the autotrophic and heterotrophic components of the ecosystem. Specifically, the macroalgae exhibit positive GPP as they photosynthesize, grow and increase in abundance, but this results in a concomitant increase in heterotrophic metabolism, which enhances Re. In the face of deforestation, both GPP and Re are reduced, resulting in little to no changes in NEP. Thus, we propose that, GPP, Re and the range between them are better measures of changes to primary productivity than NEP. Combining these with estimates of macroalgal and invertebrate diversity and abundance revealed that the Aleutian Archipelago suffered geographically variable losses to ecosystem function following widespread deforestation.

Materials and methods
While many past experiments examining primary production by autotrophic communities have relied on laboratory experiments that do not incorporate natural fluctuations in abiotic conditions, recent studies have identified techniques that measure primary production in situ, thereby increasing the ecological realism of their experiments [54][55][56][57]. For example, in situ chamber designs have been developed for estimating primary production by individual species [55,56] and whole benthic communities [29,56,57]. In general, estimates of net ecosystem production (NEP) on the benthos can be made by measuring changes in dissolved oxygen within chambers that are placed in situ over macroalgae and invertebrate communities. In this study, we deployed collapsible benthic isolation tents (cBITs) modelled after those described by Haas et al. [58] and Calhoun et al. [59] that directly measured in situ benthic oxygen production and allowed us to estimate gross primary production (GPP), ecosystem respiration (Re) and net ecosystem production NEP by the benthic communities [28,29,55]. These cBITs were the same ones used by Sullaway and Edwards [60] to measure loss of primary productivity following the

Experimental design
Our cBITs were made from 0.106 cm polycarbonate plastic triangle sheets glued to fiberglassreinforced vinyl panels (Fig 8). The frames were reinforced using stainless steel tubes with stainless steel cable to facilitate handling and to ensure they held their pyramidal shape with an internal volume of 192 L and a basal area of (0.64 m 2 ). The cBITs each had 26" skirts around the perimeter, upon which chain was laid to hold them to the benthos and prevent water exchange with the surrounding environment. This was verified by injecting fluorescein dye into the chambers and examining the perimeters for leaks. The polycarbonate walls were thin and flexible to allow hydrodynamic energy transfer into the cBITs, thereby reducing boundary layer formation around the macroalgal thalli. We verified this energy transfer using dissolving plaster blocks placed within cBITs, and by using video analysis of internal seaweed and fluorescein dye movements within the chambers relative to seaweeds outside them [60]. Sensor arrays that included a Photosynthetic Active Radiation (PAR) sensor (Odyssey Dataflow Systems Ltd), and a Dissolved Oxygen (DO mg/L) and Temperature (˚C) sensor (MiniDOT Logger, PME) were placed at the center of each cBIT (Fig 8).
During two cruises aboard the R/V Oceanus in 2016 and 2017, we deployed cBITs in each of the three habitats (kelp forest, urchin barrens, transition habitats) on each of nine islands (Table 1; Figs 1,2 and 8) for 36-hour periods to measure both day and night patterns of NEP and Re, and to ensure we captured a complete diurnal cycle. These islands span more than 1000 km and therefore experience differences in temperature, salinity, wave exposure and other biotic factors [62]. Consequently, all cBITs deployments were done in the summer (i.e. July) of each year, in similar depths (i.e. 6-8 m), and under similar wave exposures (i.e. protected from ocean swells) in order to standardize factors that could affect productivity measurements. The three habitat types were selected based on non-overlapping community assemblages (i.e., kelp forests were chosen based on abundant E. fistulosa and dense assemblages of understory macroalgae; transition habitats were chosen based on abundant E. fistulosa, little-to-no understory macroalgae, and high abundances of urchins; urchin barrens were chosen based on no E. fistulosa, little-to-no understory macroalgae and abundant urchins). These were then grouped in each island to reduce the effects within-island spatial heterogeneity in other environmental factors. For each deployment, three replicate cBITs were placed on the benthos over targeted assemblages within each habitat type. However, occasionally, replicates were lost due to logistical difficulties associated with the chamber-benthos seals ( Table 1). The water within each cBIT was replaced once per day by opening the side of the chamber and completely replacing the water with new ambient seawater to reduce "chamber effects" (i.e., the build-up of oxygen and depletion of inorganic carbon and nutrients). After each deployment, the chambers and sensors were retrieved. At six of the islands (Table 1), all organisms within each of the chambers' benthic footprints were collected, brought back to the ship, enumerated and weighed during our 2017 cruise. We measured NEP over the whole diurnal cycle, Re during the nighttime hours, and calculated GPP during the day for each cBIT during each incubation period separately according to Olivé et al. [57]. Specifically, measurements made during the night (the dark) were used to infer rates of Re, which were then combined with measurements of NEP to estimate GPP by the autotrophs [18][19][20]. Ethical Approval: All procedures performed in studies involving fishes were in accordance with the ethical standards of the institution or practice at which the studies were conducted (University of Alaska Fairbanks Institutional Animal Care and Use Committee; Permit Number: 899401-4).

Statistical analyses
All analyses were done in either Systat ver. 12 or Primer ver 6. Prior to analyses, all data were evaluated for normality by graphical examination of the residuals, which suggested they were slightly non-normal. Data for NPP, GPP, Re and the Range between GPP and Re were then square-root transformed and re-graphed, which suggested the problems were corrected. Data for PAR were log transformed, which corrected the problem. The transformed data were then examined for equality of variances using Bartlette's tests, which indicated they were homoscedastic. We then evaluated if urchin biomass, PAR, GPP, Re, NEP and the range between GPP and Re varied among the three habitats (kelp forests, urchin barrens, and transition habitats), the nine islands, and between the two study years using separate three-way Model III Nested ANOVAs, with year and habitat type as fixed factors, island nested within year as a random factor. We evaluated if the relationship between GPP and Re varied among habitats using ANCOVA, with Re as the response variable, GPP as the covariate, and habitat type as the categorical independent variable. We evaluated if the ratios in any of the habitats differed from 1.0 (i.e. GPP = Re) by assessing if the value 1.0 occurred within the 95% confidence intervals around their average values.