Stable isotopes track the ecological and biogeochemical legacy of mass mangrove forest dieback in the Gulf of Carpentaria, Australia

A combination of elemental analysis, bulk stable isotope analysis (bulk SIA) and compound-specific stable isotope analysis of amino acids (CSIA-AA) was used to assess and monitor carbon (C), nitrogen (N) and sulfur (S) cycling of a mangrove ecosystem that suffered mass dieback of trees in the Gulf of Carpentaria, Australia in 2015–2016, attributed to an extreme drought event. Three field campaigns were conducted 8, 20 and 32 months after the event over a period from 2016 to 2018 to obtain biological time-series data. Invertebrates and associated organic matter including mangroves and sediments from the impacted ecosystem showed enrichment in ¹³C, ¹⁵N and ³⁴S relative to those from an adjacent unimpacted reference ecosystem, likely indicating lower mangrove carbon fixation, lower nitrogen fixation and lower sulfate reduction in the impacted ecosystem. For example, invertebrates representing the feeding types of grazing, leaf feeding and algae feeding were more ¹³C enriched at the impacted site, by 1.7 ‰–4.1 ‰, and these differences did not change over the period from 2016 to 2018. The CSIA-AA data indicated widespread ¹³C enrichment across five essential amino acids and all groups sampled (except filter feeders) within the impacted site. The seedling density increased from 0.2 m⁻² in 2016 to 7.1 m⁻² in 2018 in the impacted forest, suggesting recovery of the vegetation. Recovery of CNS cycling, however, was not evident even after 32 months, suggesting a biogeochemical legacy of the mortality event. Continued monitoring of the post-dieback forest is required to predict the long-term trajectory of ecosystem recovery. This study shows that time-series SIA can track biogeochemical changes over time and evaluate recovery of an impacted ecosystem from an extreme event.

1 Introduction

Low-frequency, high-intensity weather events, such as droughts, tropical cyclones, heatwaves and climatic extremes can cause mass mortality of foundation species such as mangroves (Sippo et al., 2018), saltmarshes (Silliman et al., 2005), seagrasses (Thomson et al., 2015), kelps (Wernberg et al., 2016) and corals (Hughes et al., 2017). The frequency and intensity of extreme climatic events are expected to increase due to climate change (Coumou and Rahmstorf, 2012; Stott, 2016). In 2015–2016, an extensive area (>7000 ha) of mangrove forest along ∼1000 km of coastline in the Gulf of Carpentaria, Australia, experienced severe dieback, an event associated with the climatic extreme of drought (Duke et al., 2017; Sippo et al., 2020a). As mangroves show characteristics of pioneer species (Tomlinson, 2016), large-scale disturbances have likely played an important role in their evolution. However, the processes, rates and patterns of recovery from disturbances are still largely unknown. In most cases, recovery of mangroves primarily relies on the recruitment of seedlings (Smith et al., 1994; Krauss and Osland, 2019). Disturbances in mangrove forests not only affect recruitment but can also change the cycling of C, N and S. Loss of mangrove trees and root structures can change organic matter inputs, sediment oxygenation and degradation of sediment organic matter, with consequences for benthic assemblages (Sweetman et al., 2010; Bernardino et al., 2018; Harada et al., 2019), the coastal carbon cycle (Jeffrey et al., 2019; Sippo et al., 2020b), sediment C and N stocks (Adame et al., 2018), microbial assemblages, and associated nutrient processes, e.g. nitrogen fixation and sulfate reduction (Sjöling et al., 2005). Although mangroves can recover from mortality events, the rate of recovery can be slow. For example, a study of mangrove mortality attributed to an oil spill incident shows that full canopy recovery may take over 50 years (Connolly et al., 2020). Full recovery of below-ground C and N stocks may take over 40 years after mangrove replantation (Adame et al., 2018).

Stable isotope analysis (SIA) provides biogeochemical process information integrated over time and is useful for environmental assessment and monitoring. As elements such as carbon (C), nitrogen (N) and sulfur (S) circulate in the biosphere, stable isotopic compositions of ¹³C ∕ ¹²C, ¹⁵N ∕ ¹⁴N and ³⁴S ∕ ³²S can change in predictable ways due to mixing and fractionation, giving insights into sources and cycling of these elements (Fry, 2006). SIA has been widely used in mangrove ecosystem studies to better understand food web interactions (Bouillon et al., 2008; Larsen et al., 2012; Bui and Lee, 2014; Abrantes et al., 2015); nutrient uptake (McKee et al., 2002); water use (Santini et al., 2015; Hayes et al., 2019); cycling of C (Maher et al., 2013a, 2017; Sasmito et al., 2020), N (Fry and Cormier, 2011) and S (Raven et al. 2019); and greenhouse gas emissions (Maher et al., 2013b). While traditional field methods such as measuring species composition to evaluate the structure and functioning of ecosystems can be time-consuming and expensive, SIA of ecosystem components can evaluate functional aspects of element cycling and food webs in a cost-effective way (Fry, 2006). To quantify food web dynamics, SIA of total organic matter (“bulk”) requires determination of the baseline isotope values of the food web, but this is difficult to achieve, particularly for complex detrital food webs. Some of the uncertainties associated with bulk SIA have been clarified by compound-specific isotope analysis of amino acids (CSIA-AA). This technique has been increasingly used to assess pathways of energy transfer throughout food webs by distinguishing the effects of baseline isotope values and trophic transfer. Distinct isotopic fractionation between two groups of amino acids occurs with each trophic transfer. In general, non-essential and/or trophic AAs show large isotopic fractionation per trophic step, while essential and/or source AAs show little fractionation, reflecting the baseline isotope values of the food web (Ishikawa et al., 2018; Larsen et al., 2013; Ohkouchi et al., 2017).

We investigated changes in C, N and S cycling associated with the Gulf of Carpentaria mangrove forest dieback (Duke et al., 2017), using a combination of traditional ecological survey techniques, bulk SIA and CSIA-AA of carbon. We hypothesised that the mortality of mangrove foundation species has changed the overall circulation of C, N and S elements and that these biogeochemical changes would most likely be reflected in δ¹³C, δ¹⁵N and δ³⁴S values of mangrove ecosystem components such as mangrove plants, sediment and associated animals. We also tested the hypothesis that these isotopic compositions changed over time with the recovery of mangrove vegetation. δ¹³C, δ¹⁵N and δ³⁴S values were measured for samples including mangroves, sediment and invertebrates collected in a comparative setting of an impacted mangrove forest site and an adjoining unaffected reference forest site in the Gulf of Carpentaria, Australia.

2 Material and methods

2.1 Study site

The Gulf of Carpentaria in tropical northern Australia is an extensive, shallow coastal gulf. The area mainly consists of low-lying wetlands and is largely inaccessible with little direct human activity. Mangroves are abundant in the area, but the dry climate limits the extent, diversity and height of mangroves in the region (Asbridge et al., 2016). The wide tidal wetlands spread along the shoreline with high intertidal saltpans and saltmarsh covering more area than mangroves. Mangroves are often distributed in the seaward margin, typically as a narrow strip and fronted by extensive, shallow mudflats. The distribution of mangroves in this region is associated with tidal and freshwater inundation, river discharge and regular sediment supply through freshwater input. Increased amounts of rainfall and associated flooding and sea level rise were responsible for recent mangrove extension in this region between 1987 and 2014 (Asbridge et al., 2016). Avicennia marina was the dominant mangrove species at the study site in Karumba, Gulf of Carpentaria (Fig. 1a).

Figure 1 The study location at Karumba in the Gulf of Carpentaria, Queensland, Australia (−17.435572^∘ S, 140.844766^∘ E; image provided by ©Google Earth). (a) Three sampling transects within the unimpacted reference site and three within the impacted site (shown as a white line). (b, c) Representative transects from the unimpacted (b) and impacted (c) sites. Each transect was approximately 200 m. Samples were collected along each transect from higher to lower intertidal zones. Red indicated mass dieback region.

Over 7000 ha of mangroves along ∼1000 km of the Gulf of Carpentaria coastline in Australia experienced mass mortality during the summer in 2015–2016 (Duke et al., 2017), the most extensive mangrove forest dieback ever recorded due to natural causes (Sippo et al., 2018). At the same time, there were coinciding mangrove mass mortality events in Exmouth, Western Australia (Lovelock et al., 2017), and Kakadu National Park, Northern Territory (Asbridge et al., 2019). The climate in the Gulf region is wet–dry tropical with mean annual precipitation ranging from approximately 600 to 900 mm. Dry conditions prevail for 6 to 8 months, and most rainfall occurs between December and March (Bureau of Meteorology; see http://www.bom.gov.au, last access: November 2019). The dieback coincided with a weak monsoon (low rainfall), combined with a high vapour pressure deficit and El Niño–Southern Oscillation-induced low sea levels (Duke et al., 2017; Lovelock et al., 2017; Harris et al., 2017). The drought conditions most likely caused accumulative hydric, thermal and radiant stresses (Duke et al., 2017). In addition to low water availability, iron (Fe) toxicity due to a rapid mobilisation of sedimentary Fe and regional variability in groundwater flows may have also played a role in the dieback (Sippo et al., 2020a). The event led to the widespread death of mangrove trees in the region, providing an unfortunate yet unique opportunity to test tree mortality effects on the biogeochemical and ecological functioning of mangroves and capture recovery patterns.

Three field campaigns were conducted in August 2016 (8 months after the event), August 2017 (20 months after the event) and August 2018 (32 months after the event) in the winter dry seasons in Karumba. Some local factors (e.g. river influence and localised groundwater flow paths) may have kept some mangroves from dying back in the region. A forest that had suffered dieback (impacted) on the east side of the Norman River outlet and an adjoining unaffected forest (unimpacted) on the west side provided the setting for comparisons. In order to assess differences between the two forests (impacted vs. unimpacted), as well as to capture trends from across the intertidal zone and to ensure that the physical–oceanographic conditions between the two forests were as similar as possible, three sampling transects (2 to 2.5 km apart) were set for each forest, with the length of each transect being approximately 200 m (Fig. 1b, c). Each transect consists of six sampling plots (approx. 40 m apart), namely forest edge (landward), high, mid, low, forest edge (seaward) and mudflat (Table 1). Samples from each plot were pooled for analysis. Due to logistical constraints and the presence of saltwater crocodiles, fieldwork was restricted to daytime, low tide and dry seasons.

Table 1 Spatial and temporal details of the sampling design (x = sampled).

^* For invertebrate SIA, to gain sufficient sampling size, high, mid and low plots were pooled for analysis.

2.2 Samples

Since our focus was to measure the recovery of mangrove vegetation and food webs, we monitored mangrove saplings and seedlings and stable isotopes of invertebrates during the period from 2016 to 2018. Mangrove and sediment samples were also collected, but they are limited to 2018 (Table 1). Some of the SIA samples including invertebrates, mangrove and sediment collected in 2017 were used to measure the initial dieback reported in Harada et al. (2019). During each field campaign, four common mangrove macroinvertebrate groups with different feeding modes were collected from each forest including a leaf-eating crab (Parasesarma moluccensis and/or Episesarma sp.), an algae-eating (deposit-feeding) crab (Tubuca signata), a grazer gastropod (Telescopium telescopium) and a filter-feeding bivalve (Saccostrea sp., an oyster). For each feeding group, three to five individuals at each of the sampling transects (n=3) within the forest were collected and muscle tissues were pooled for SIA. In 2018, fully developed green leaves of A. marina were collected from about 1 to 1.5 m height, from three to five individual trees (one to three leaves per tree) at each sampling plot, stored in plastic containers and then composited. In the impacted site, leaves were collected from regrowth of trees that had survived. Leaf samples were washed thoroughly and rinsed with distilled water, and the main vein was removed. Additionally, wood samples (n=2, 5 to 25 cm diameter) were collected using a hand saw from stems at chest height from the mid intertidal zone of each forest. Dead trees were sampled at the impacted site. Two to three bulk SIA measurements were made from the sapwood (2 to 3 cm deep) of each sample and measurements were averaged. Wood samples were generally very low in S and not sufficient for δ³⁴S analysis.

In 2018, surface (<0.5 cm) sediments that represent the recent deposition and microphytobenthos (MPB) were collected along each transect. Additionally, subsurface (0.5 to 20 cm) sediment samples (n=6) that represent long-term averages were collected at the mid intertidal zone of each forest using a core sampler (5 cm in diameter and 20 cm deep). For δ¹³C measurements, the sediment samples were acidified with 1 M HCl to remove the inorganic fraction. MPB samples (n=6) were separated from surface sediment collected at each forest. The separation was achieved through density gradient centrifugation in Ludox colloidal silica (Sigma) as described in Bui and Lee (2014). The MPB extraction was followed by microscopic examination to confirm that samples mostly contained green cells (i.e. diatoms and filamentous cyanobacteria). Additionally, surface sediment samples were collected offshore (approx. 1 km) using a grab sampler and from the adjacent saltpan (approx. 200 m inland from the forest). Offshore water samples (n=3) were also collected and then filtered through glass fibre filters (diameter 44 mm, pore size 0.7 µm, Whatman GF/F) to obtain particulate organic matter (POM), material such as phytoplankton that is transported to the mangrove.

To estimate mangrove seedling and sapling densities (individuals per square metre) from each forest and their changes over time, seedlings and saplings were counted with a 50 cm × 50 cm quadrat at the mid intertidal zone. A photo was taken of each quadrat (for 2016, n=124 for the unimpacted forest and n=143 for the impacted forest; for 2017, n=161 and n=175; and for 2018, n=80 and n=117, respectively), and then counts of seedlings and samplings were made in the laboratory. The seedlings and samplings were mostly A. marina but also include some Aegiceras corniculatum.

2.3 Stable isotope analysis

All invertebrate, leaf, sediment and filter samples were stored separately in sealed plastic containers at −20 ^∘C until analysis, then dried at 60 ^∘C, powdered, homogenised and put in tin capsules for SIA. δ¹³C, δ¹⁵N and δ³⁴S measurements were carried out on an elemental analyser (Europa EA-GSL, Sercon) coupled to an isotope ratio mass spectrometer (Hydra 20-22, Sercon) at Griffith University, Brisbane, Australia. Isotope values are reported relative to Vienna Pee Dee Belemnite (PDB), atmospheric N₂ (AIR) and Vienna Canyon Diablo Troilite (VCDT) for C, N and S, respectively. Harada et al. (2019) reported δ¹³C and δ¹⁵N values for some of the samples collected in 2017 from the same study location.

The samples collected in 2017 showed substantial differences in δ¹³C values between the two forests (Harada et al., 2019), suggesting that they are representative of the dieback impact. The 2017 samples from each forest including the mangrove leaf (n=2), MPB (n=1), algae feeder (n=3), leaf feeder (n=2 to 3), grazer (n=3), and filter feeder (n=2) were further measured for carbon isotopic composition of individual amino acids (δ¹³C_AA). A total of 8 mg (for animal tissues) or 30 mg (for plant tissues) of sample materials were transferred to borosilicate vials with heat- and acid-resistant caps. They were then flushed with N₂ gas, sealed and hydrolysed in 0.5 mL (animal tissues) or 2 mL (plant tissues) of 6 M HCl at 150 ^∘C for 70 min, then dried in a heating block at 60 ^∘C under a stream of N₂ gas. The dried samples were derivatised by methoxycarbonylation as described by Walsh et al. (2014). Amino acid derivatives were separated by a Trace GC Ultra gas chromatograph (GC) (Thermo Scientific) using a DB-23 column (Agilent, 30 m × 0.25 mm, 0.25 µm film) at the stable isotope facility at the University of California (Davis, CA, USA). The GC was interfaced with a Delta V Plus isotope ratio mass spectrometer via a GC IsoLink (Thermo Scientific). L-Norleucine was used as an internal standard and to calculate provisional values. Pure AA mixtures with calibrated δ¹³C were co-measured. One mixture was used for final calibration and others were for the scale normalisation standard and the primary quality assurance standard (unused in corrections). Two working standards were co-measured as secondary quality assurance materials. Exogenous carbon was accounted for using the method detailed by Docherty et al. (2001). Following these processes, δ¹³C values were determined for 10 AAs (Gly, glycine; Asx, aspartic acid/asparagine; Pro, proline; Glx, glutamic acid/glutamine; Ala, alanine; Lys, lysine; Ile, isoleucine; Leu, leucine; Phe, phenylalanine; and Val, valine). Met, methionine; His, histidine; and Hyp, hydroxyproline were at or below the limit of quantitation (LOQ) for some samples. Since we were interested in δ¹³C values of essential amino acids (δ¹³C_EAA), δ¹³C values of Lys, Ile, Leu, Phe and Val were reported in the study.

2.4 Data analysis

All statistical analyses were undertaken in R version 3.4.3 with RStudio interface version 1.1.414. Differences among group means were explored with ANOVA, but for the count data, i.e. seedling and sampling density, a generalised linear model (GLM) with Poisson distribution was used. Before performing ANOVA, the assumptions of homogeneity of variance and normality were tested using Levene and Shapiro–Wilk tests, respectively. Two-way ANOVA was used to test the effects of time (year) and forest type (unimpacted and impacted) on stable isotope values of invertebrates. To explore δ¹³C patterns among five EAAs, δ¹³C_EAA values were normalised to the respective sample means following the procedure of Larsen et al. (2009) as follows:

$$\begin{array}{}\text{(1)}& \mathrm{Norm}\left({\mathit{\delta }}_{\mathrm{EAA}}\right)={\mathit{\delta }}_{\mathrm{EAA}}-\mathit{\mu },\end{array}$$

where μ represents the mean value of all five EAAs (Ile, Leu, Lys, Phe and Val) in the sample. PERMANOVA was performed to test if the pattern of δ¹³C among five EAAs of samples differ between the forests. In this analysis, the normalised δ¹³C_EAA dataset was used and the Euclidean distance was used as distance metric. Permutation test of the multivariate homogeneity of dispersions was performed to check whether dispersions around the centroids are similar between the two forests. All statistical tests used a significance criterion of α=0.05.

3 Results

3.1 Mangroves

In 2016, at the impacted site, mangrove seedling and sapling populations were lower than at the unimpacted site (GLM df  = 1, estimate 3.75, p<0.001) but significantly increased throughout the period from 2016 to 2018 (GLM df = 2, estimate 0.01, p<0.001; Fig. 2). C, N and S elemental compositions (%) of A. marina did not differ significantly between the two sites, but the isotopic compositions varied significantly (Table 2). The δ¹³C values (mean ± SD) of green leaves harvested from A. marina trees were significantly higher (i.e. more ¹³C enriched) in the impacted site ($-\mathrm{25.8}\pm \mathrm{1.0}$ ‰) than the unimpacted site ($-\mathrm{28.4}\pm \mathrm{1.5}$ ‰; ANOVA $F_{\mathrm{1},\mathrm{28}}=\mathrm{32.9}$, p<0.001; Fig. 3a, Table 2). The δ¹⁵N values varied more in the impacted site (ranged from −0.9 ‰ to 6.7 ‰) than in the unimpacted site (ranged from 2.9 ‰ to 6.2 ‰; Fig. 3a). The δ³⁴S values were generally higher in the impacted site (13.5±5.4 ‰, range 7.7 ‰ to 23.3 ‰) than the unimpacted site (12.6±5.6 ‰, range 5.0 ‰ to 21.9 ‰). Leaf δ³⁴S values became more ³⁴S depleted from higher to lower intertidal zones in the impacted site (ANOVA $F_{\mathrm{4},\mathrm{10}}=\mathrm{5.56}$, p=0.013; Fig. 3b). This pattern was weaker in the unimpacted site, but leaf δ³⁴S values significantly varied across intertidal zones (ANOVA $F_{\mathrm{4},\mathrm{10}}=\mathrm{6.48}$, p=0.007; Fig. 3b). Leaf δ¹³C and δ¹⁵N values did not display such patterns along the intertidal zones (Table S2 in the Supplement). Leaf C, N and S (%) did not show any clear trends among the two forests and along transects. Yellow leaves generally had a higher S content (∼1.2 %) than green leaves (0.5 % to 0.7 %) (Table 2). Wood samples generally had very low N and S contents, significantly lower than the leaves.

Figure 2 Recovery of mangrove vegetation at the impacted site during a 2-year period from 2016 to 2018 (a and b; approx. 8 and 32 months, respectively, after the dieback event). Seedling and sampling densities of mangrove species (mostly, A. marina) significantly increased in the impacted site (c; Table S1).

Figure 3 CNS isotopic compositions of green leaves of A. marina from the unimpacted and impacted sites. All samples were collected in 2018, 32 months after the dieback. (a) Leaf δ¹³C and δ¹⁵N values. (b) Leaf δ³⁴S values across the intertidal zones. Error values are SE.

Table 2 Elemental and isotopic compositions of the mangrove A. marina (mean, SD) in 2018, 32 months after the dieback.

3.2 Sediment

For the surface (< 0.5 cm) sediment collected in 2018, TOC (%) differed significantly between the two forests, with the values (mean ± SD) of 2.02 ± 1.16 % for the unimpacted site and 1.06 ± 0.37 % for the impacted site (ANOVA $F_{\mathrm{1},\mathrm{28}}=\mathrm{12.75}$, p=0.001), suggesting that the surface sediment from the impacted site contains ∼ 48 % lower TOC relative to that of the unimpacted site (Table 3). The surface sediment TN (%) was also significantly lower for the impacted site (0.09 ± 0.03) than the unimpacted site (0.15 ± 0.06 %; ANOVA $F_{\mathrm{1},\mathrm{28}}=\mathrm{9.32}$, p=0.005). TOC of mudflat (< 0.5 cm) sediment collected adjacent to the two forests also differed significantly, with those of the impacted site being lower (0.58 ± 0.18 %) than the unimpacted site (1.02 ± 0.08 %; ANOVA $F_{\mathrm{1},\mathrm{4}}=\mathrm{14.54}$, p=0.019). The TN (%) of mudflat (< 0.5 cm) sediment was also significantly lower for the impacted site (ANOVA $F_{\mathrm{1},\mathrm{4}}=\mathrm{9.81}$, p=0.035). The TOC (%) of the 0.5–20 cm sediment did not differ significantly between the two sites, with the values of 1.83 ± 0.73 % for the unimpacted and 1.29 ± 0.55 % for the impacted site (ANOVA $F_{\mathrm{1},\mathrm{10}}=\mathrm{2.07}$, p=0.181). The TN (%) of the 0.5–20 cm sediment also did not differ significantly between the two sites, with the values of 0.09 ± 0.03 % for the unimpacted and 0.12 ± 0.04 % for the impacted site (ANOVA $F_{\mathrm{1},\mathrm{10}}=\mathrm{3.02}$, p=0.113). Despite the substantial variation in TOC and TN, the C∕N ratio did not differ significantly between the two sites (ANOVA p>0.05).

Table 3 Elemental and isotopic compositions of sediment in 2018, 32 months after the dieback (mean, SD).

^* Values were taken from Harada et al. (2019).

The δ¹³C values of surface sediment differed significantly between the two sites, with those from the impacted site showing higher values (−21.8 ± 1.0 ‰) than those from the unimpacted site (−24.3 ± 1.2 ‰; ANOVA $F_{\mathrm{1},\mathrm{28}}=\mathrm{22.48}$, p < 0.001). This pattern was consistent across the intertidal zones, with δ¹³C values from the impacted site becoming similar to those from the adjacent mudflat (Fig. 4; Table S3). However, those of the 0.5 to 20 cm sediment did not differ significantly, with the values of −24.4 ± 0.5 ‰ for the impacted site and −25.2 ± 0.9 ‰ for the unimpacted site (ANOVA $F_{\mathrm{1},\mathrm{10}}=\mathrm{3.92}$, p=0.076). Surface sediment collected in the adjacent mudflat did not display a significant difference in δ¹³C values between the two sites, with the values of −21.2 ± 0.9 ‰ for those collected adjacent to the unimpacted forest and −21.8 ± 0.8 ‰ for those collected adjacent to the impacted site (ANOVA $F_{\mathrm{1},\mathrm{4}}=\mathrm{0.64}$, p=0.47). The δ¹³C value of the surface sediment (−21.8 ± 1.0 ‰) in the impacted forest was similar to the value of those collected in the mudflat (−21.2 ± 0.9 ‰), which were also similar to those collected from offshore (−21.5 ± 1.1 ‰). Those δ¹³C values also matched with the δ¹³C value of POM collected offshore (−21.5 ± 1.5 ‰). Surface sediment collected from adjacent unvegetated saltpan areas also showed similar values (Table 3). MPB extracted from the surface sediment showed significantly different δ¹³C values between the impacted site (−21.5 ± 1.3 ‰) and unimpacted site (−25.2 ± 1.0 ‰; ANOVA $F_{\mathrm{1},\mathrm{10}}=\mathrm{28.53}$, p < 0.001).

Figure 4 C elemental and isotopic compositions of surface (< 0.5 cm) sediment along the unimpacted reference transects vs. impacted transects (n=3 per data point). All samples were collected in 2018, 32 months after the dieback. Error values are SE. (a) Sediment TOC, %. (b) Sediment δ¹³C values, ‰.

3.3 Fauna

CNS isotopic compositions of mangrove macroinvertebrates representing algae feeder, grazers and leaf feeders from the impacted site were more enriched in ¹³C, ¹⁵N and ³⁴S than their counterparts from the unimpacted forest throughout the period between 2016 and 2018 (Fig. 5). However, the filter-feeding oyster that relies on water column organic matter showed relatively less difference between the two forests (Fig. 5). Overall, δ¹³C values of the four feeding groups range from −23 ‰ to −15 ‰ for the unimpacted site and −20 ‰ to −14 ‰ for the impacted site. The δ¹⁵N values ranged from 5.5 ‰ to 9.1 ‰ for the unimpacted site, with the fauna in the impacted site having a slightly higher range of 5.6 ‰ to 9.5 ‰. The δ³⁴S values ranged from 8.2 ‰ to 16 ‰ for the unimpacted site and 13.4 ‰ to 21.7 ‰ for the impacted site. The δ¹³C, δ¹⁵N and δ³⁴S values of the algae feeder and the grazer significantly differed between the two forests (ANOVA p < 0.05). The δ¹³C and δ³⁴S values of the leaf feeder significantly differed between the two forests (ANOVA p < 0.05), but the δ¹⁵N values were not significantly different (ANOVA $F_{\mathrm{1},\mathrm{13}}=\mathrm{1.72}$, p = 0.212). The filter feeder δ¹³C values did not differ between the two forests (ANOVA $F_{\mathrm{1},\mathrm{8}}=\mathrm{1.719}$, p=0.212), but the δ¹⁵N and δ³⁴S values differed significantly between the two forests (ANOVA p < 0.05). The effect of time was significant for the leaf feeder δ³⁴S values, the grazer δ¹⁵N values, the filter feeder δ¹⁵N and δ¹³C values (ANOVA p < 0.05). In 2018, leaf feeder δ¹³C values, grazer δ³⁴S values and algae feeder δ³⁴S values significantly differed between the impacted and unimpacted sites (Turkey post hoc test, p < 0.05), showing no recovery of the invertebrate fauna δ¹³C and δ³⁴S status 32 months after the dieback event. However, δ¹⁵N values became similar between the two forests in 2018 (Fig. 5).

Figure 5 Changes in CNS isotopic compositions of mangrove macrofaunal groups with four different feeding modes from 2016, 2017 and 2018 (i.e. 8, 20 and 32 months after the event) between the unimpacted reference and impacted mangrove forest sites. Error bars are ± SE (n=2 to 6 per data point). An asterisk (^*) indicates a significant difference in the year. Mean ± SE values and sample sizes are also provided in Table S4.

3.4 Compound-specific isotope analysis of amino acid carbon

The samples collected in 2017 were further measured for carbon isotopic compositions in individual essential amino acids. δ¹³C_EAA values corresponded to the bulk δ¹³C values, with the samples from the impacted site consistently showing higher values than those from the unimpacted forest (Table 4 and Fig. 6). The pattern of normalised δ¹³C among five EAAs (Lys, Ile, Val, Leu and Phe) for all the consumers did not differ between the two forests (p>0.05; Table 4 and Fig. S1 in the Supplement). The δ¹³C_EAA pattern of mangrove leaves did not differ between the forests, regardless of the substantial bulk isotope difference between the unimpacted (−26.7 ± 2.2 ‰) and impacted (−25.4 ± 0.1 ‰) forests (Table 4). This isotope pattern was similar for all four feeding groups as well as for the MPB samples (Figs. 6 and S1).

Figure 6 C isotopic compositions in essential amino acids (EAAs) for four mangrove consumer groups and resources including mangrove leaves (A. marina) and MPB from the unimpacted and impacted mangrove sites during 2017 (20 months after the dieback). While there are clear offsets in individual δ¹³C_EAA values between the two forests, normalised δ¹³C_EAA fingerprint patterns as per Larsen (2009) shown in Fig. S1 did not differ (PERMANOVA p>0.05, Table 3). Error bars show ± SD. The data are provided in Table S5.

Table 4 Bulk δ¹³C values, mean δ¹³C values of five EAAs (‰ ) and differences (Δ, ‰) between the two forests in 2017, 20 months after the dieback (mean ± SD).

4 Discussion

4.1 Mangroves

The recovery of mangrove forests from tree mortality events caused by disturbances such as cyclones generally relies upon the recruitment of seedlings (Smith et al., 1994; Krauss and Osland, 2019). Subsequently, degraded habitats with a reduced seed pool, production and delivery, e.g. by habitat fragmentation, may show slower forest recovery (Milbrandt et al., 2006). The establishment of seedlings may also be inhibited by persistent inundation due to a decreased sediment elevation (Cahoon et al., 2003; Asbridge et al., 2018). Mangroves are resilient ecosystems and may recover quickly from natural disturbances (Sherman et al., 2001), but in some cases, full recovery may take decades (Imbert et al., 2000; Connolly et al., 2020). Further, mangrove degradation may be followed by fast colonisation of non-mangrove herbaceous species, e.g. succulent saltmarsh (Mbense et al., 2016; McKee et al., 2007; Rashid et al., 2009). However, this was not largely evident from our study site. In both mangrove forests at the Gulf of Carpentaria site, the density of mangrove seedlings and samplings significantly increased throughout the period from 2016 to 2018, suggesting that recovery was starting to occur in some areas within 32 months after the dieback and propagule pool was available in the vicinity. The increase in seedling and sampling density at the unimpacted site was unexpected, but this indicates that there was some stress at the unimpacted site during the dieback period, and/or the temporal variability of seedling and sampling density was high at the site.

The substantial differences in CNS isotopic compositions in A. marina occurring between the two sites suggested differences in the environmental conditions and biogeochemical processes that were possibly associated with the mangrove mortality effect. The leaf δ¹³C values in the impacted forest were relatively enriched in ¹³C. This C isotope pattern may be due to reduced stomatal conductance that causes lower internal carbon dioxide concentrations and lower carbon isotope fractionation (Farquhar et al., 1989; Lin and Sternberg, 1992a; Lin and Sternberg, 1992b). Higher leaf δ¹³C values can also be associated with increased carboxylation efficiency associated with higher nutrients, e.g. N in leaves (Cordell et al., 1999) and thicker leaves with higher internal resistance to carbon dioxide diffusion. Younger leaves can show higher δ¹³C values than aged leaves due to ¹³C enriched fractions (e.g. carbohydrates) transported from older autotrophic leaves to more heterotrophic young leaves (Werth et al., 2015). Leaves exposed to full sun can show higher δ¹³C values than shaded leaves (Farquhar et al., 1989). δ¹³C values can also vary among different plant tissues in A. marina (Kelleway et al., 2018). Leaf N (%) did not differ between the two forests, suggesting that the two sites may have similar plant N availability. Previous studies show additions of nutrients such as N and/or P did not play a considerable role in mangrove leaf δ¹³C variations (McKee et al., 2002), but salinity played an important role (Lin and Sternberg, 1992b). For example, leaf δ¹³C values of A. marina at a lower salinity site were relatively depleted in ¹³C and averaged about −31 ‰ (Kelleway et al., 2018). Overall, higher leaf δ¹³C values in the impacted forest likely suggest that there are chronic stresses associated with the dieback event that reduced stomatal conductance. Such environmental stresses may include hypersalinisation of sediments and hydric, thermal and radiant stresses following canopy losses that cause higher evaporation and lower water availability. The leaves at the unimpacted site were largely depleted in ¹³C, suggesting that there was higher water availability at the unimpacted site, possibly associated with regional variability in groundwater flows, e.g. Sippo et al. (2020).

Leaf δ¹⁵N values varied more in the impacted (ranged from −0.9 ‰ to 6.7 ‰) than in the unimpacted forest (ranged from 2.9 ‰ to 6.2 ‰), but the means were similar (4.3 ± 1.9 ‰ for the impacted and 4.4 ± 0.8 ‰ for the unimpacted), suggesting that two sites have similar background δ¹⁵N conditions. Generally, leaf δ¹⁵N varies due to N sources, microbial processes that enrich or deplete ¹⁵N in soil or water, and isotope fractionation during plant N uptake. Previous studies showed that in pristine mangrove forests, leaf δ¹⁵N values generally range around −2 ‰ to 3 ‰ (Fry and Smith, 2002; Smallwood et al., 2003). Such low δ¹⁵N values may reflect long-term N fixation inputs (e.g. around 0 ‰) (Fogel et al., 2008) and marine nitrate inputs (Dore et al., 2002). Much higher δ¹⁵N values (>10 ‰) may be associated with anthropogenic N inputs (Fry and Cormier, 2011). Our sites showed moderate δ¹⁵N values (about 4 ‰), suggesting that in addition to N fixation inputs and marine N inputs, there may be considerable microbial ¹⁵N enrichment in dissolved inorganic nitrogen pools of ammonium and nitrate. The higher variability in leaf δ¹⁵N in the impacted forest suggests higher variability in processes affecting the δ¹⁵N status of available N. For example, changes to sediment conditions including redox transitions and soil moisture content following the dieback may have affected microbial processes of N, whereas the unimpacted forest may have had more stable sediment conditions and N processes. Isotope fractionation during plant N uptake may also be an explanation for leaf δ¹⁵N variability (Fry et al., 2000), but such fractionation is poorly known for mangroves. A study reported that additions of P nutrients increased N demand and decreased ¹⁵N fractionation (McKee et al., 2002); however, as we did not measure P, we could not determine whether this was the case.

Leaf δ³⁴S values differed considerably between the two forests, with the impacted forest generally having higher values (13.5 ± 5.4 ‰, range 7.7 ‰ to 23.3 ‰) than the unimpacted forest (12.6 ± 5.6 ‰, range 5.0 ‰ to 21.9 ‰). Leaf δ³⁴S values showed trends along the six transects, with values decreasing from the upper to lower intertidal zones. Based on previous studies, mangrove leaf δ³⁴S values generally vary between −20 ‰ and 20 ‰ (Okada and Sasaki, 1995, 1998; Fry and Smith, 2002). Higher δ³⁴S values are likely associated with seawater sulfate, which is ³⁴S enriched (i.e. 21 ‰) and due to a large isotope fractionation (up to 70 ‰) during sulfate reduction (Kaplan and Rittenberg, 1964). Lower δ³⁴S values are likely associated with sedimentary sulfide S that is ³⁴S depleted, for example, −21 ‰ (Okada and Sasaki, 1995). Leaf δ³⁴S values of around 14 ‰ to 18 ‰ suggest mangrove incorporations of seawater sulfate S (δ³⁴S, ∼ 21 ‰), with only a small isotopic fractionation occurring through absorption and assimilation steps (Okada and Sasaki, 1995). Plants generally show δ³⁴S values slightly lower than source sulfate S by an average of −1.5 ‰ (Trust and Fry, 1992). Low leaf δ³⁴S values, for instance the lowest value of 5 ‰ found in the unimpacted site, suggest that the most probable source of this ³⁴S-depleted S is sulfide oxidation, followed by mixing with seawater sulfate. Low δ³⁴S values in mangrove root vascular tissues may indicate assimilation and oxidation of sulfide, potentially to reduce their toxic sulfide exposure (Fry et al., 1982; Raven et al., 2019), with a reported isotope effect of −5.2 ‰ for the non-biological oxidation of sulfide (Fry et al., 1988) and a smaller +1–3 ‰ effect for anaerobic oxidation of sulfide by photosynthetic bacteria (Fry et al., 1984). An explanation for our observed δ³⁴S pattern may be lower plant incorporation of sulfide S in the impacted site and also in the higher intertidal zones where we expect that mangrove sediment is relatively more oxidised, and the production of sulfide may be lower due to lower sulfate reduction.

4.2 Sediment

In healthy mangrove forests, the fate of C fixed by primary producers includes burial within the sediment, atmospheric emissions and outwelling to the ocean (Maher et al., 2018), but how mangrove mortality affects such processes is poorly understood. In most cases, C within mangrove sediment decreases following forest loss due to degradation with increased CO₂ emissions (Otero et al., 2017; Adame et al., 2018). Lower TOC (%) and higher sediment δ¹³C values in the impacted forest are probably related to sediment C loss and lower autochthonous C inputs (i.e. leaf litter) following the mangrove mortality event. Consistent with this, the sediment N (%) and δ¹⁵N data showed a similar pattern, suggesting N loss and degradation. The surface sediment varied more than the subsurface fraction. This is probably because the surface fraction is generally more aerobic and remineralisation of organic matter occurs more rapidly (Burdidge, 2011). Sediment δ¹³C and ¹⁵N values can increase during degradation of sediment organic matter following mangrove loss (Adame and Fry, 2016; Adame et al., 2018). Changes in sediment C and N may also be associated with root turnover. The MPB δ¹³C values significantly differed, with those from the impacted being higher than the unimpacted. The higher values probably indicate lower respiratory inputs of CO₂ from mangroves (Maher et al., 2013b). Our findings here are consistent with the finding of Sippo et al. (2019) that changes to oceanic carbon outwelling rates following mangrove loss are likely associated with a gradual loss of sediment carbon; similar to our finding of increased sediment δ¹³C values in the impacted site, an isotope effect may have been due to loss of sediment mangrove C and/or replacement of mangrove peats with marine sediment.

4.3 Fauna

CNS isotopic compositions of consumers including an algae feeder, a grazer and a leaf feeder from the impacted site were more enriched in ¹³C, ¹⁵N and ³⁴S. These differences remained consistent throughout the three sampling periods of 2016, 2017 and 2018. Consistent with the findings from mangrove leaves, MPB and soil, the data suggested substantial changes in cycling of CNS associated with the mangrove mortality event. At the impacted site, consumers were more ¹³C enriched, likely due to the loss of ¹³C-depleted mangrove organic matter. Consumer δ¹³C values can change due to changes to available organic matter, altered feeding dependencies and changes to organic matter δ¹³C values. In our case, MPB δ¹³C values changed significantly, likely due to changes to organic matter respiratory inputs and/or altered light environment and soil moisture contents that may change isotopic fractionation during carbon fixation. The consumer δ¹³C values and their ranges at our study site are fairly consistent with the reported mangrove consumer δ¹³C values elsewhere (e.g. Lee, 2000; Bouillon et al., 2002; Demopoulos et al., 2007). The typical trophic enrichment factor for carbon isotopes in small invertebrates is about +1 ‰ (Vander Zanden and Rasmussen, 2001; McCutchan et al., 2003). Lower consumer δ¹³C values are generally associated with mangrove detritus that is depleted in ¹³C. Typical mangrove leaf-eating sesarmid crab species generally have tissue δ¹³C values within about +5 ‰ from mangrove detritus (Bui and Lee, 2014). Higher δ¹³C values of consumers are generally tied to MPB. Our MPB endmember δ¹³C values of −25.2 ‰ for the unimpacted site and −21.5 ‰ for the impacted site did not match with the consumer δ¹³C values (around −15 ‰ to −14 ‰), suggesting our characterisation of MPB endmember δ¹³C values was incomplete. This is probably because MPB can vary substantially within mangrove ecosystems (Bouillon et al., 2008), and consumers may be preferentially assimilating more ¹³C enriched fractions of MPB, for example, diatom and/or filamentous cyanobacteria that can range about −15 ‰ to −20 ‰ (Craig, 1953; Fry and Wainright, 1991). The leaf feeders were relatively depleted in ¹³C, with δ¹³C values of about −21 ‰ to −18 ‰, likely due to some use of mangrove leaves. Due to difficulties obtaining representative endmembers, mixing analysis was not achieved using these data to quantify feeding dependencies. However, MPB probably played an important dietary role in both forests because the difference in MPB δ¹³C values between the two forests were reflected in the difference in the consumers δ¹³C values between the two forests (Harada et al. 2019).

Consistent with the mangrove leave δ³⁴S values, the consumer δ³⁴S values also indicated possible changes to S cycling. The consumer δ³⁴S values were generally higher in the impacted site (range 13.4 ‰ to 21.7 ‰) than in the unimpacted site (range 8.2 ‰ to 16 ‰), suggesting lower sulfate reduction with decreased sulfide inputs at the impacted site. Fixation of sulfate by phytoplankton occurs with a small isotope effect, around 1 ‰ to 2 ‰ (Fry, 2006); therefore phytoplankton δ³⁴S values from the coastal ocean are generally close to the seawater sulfate S value of 21 ‰. MPB generally have lower δ³⁴S values than phytoplankton, e.g. near 10 ‰ for MPB in a mangrove ecosystem (Harada et al., unpublished), likely due to some use of sedimentary sulfide S (depleted in ³⁴S). Our consumer δ³⁴S values were lower than 21 ‰, suggesting some use of MPB as well as mangrove detritus. The consumer δ¹⁵N also indicates possible changes to N cycling, with the consumer in the impacted site generally having higher values than those from in the unimpacted site. The higher δ¹⁵N values are likely associated with degradation of organic matter, microbial ¹⁵N enrichment in dissolved inorganic N during degradation and lower N fixation inputs that typically show low δ¹⁵N values. While the δ¹³C and δ³⁴S values consistently differed between the two forests during the 2-year survey, the δ¹⁵N values started showing matches between the two forests in 2018, likely suggesting recovery of δ¹⁵N status to the background conditions, and/or that the recovery of N may be faster than C and S elements. This may be the case as mangrove ecosystems are generally N limited (Reef et al., 2010), and circulation of N elements is faster than those of C and S elements.

4.4 Compound-specific isotope analysis of amino acids

It is considered that environmental resources such as vascular plants and microalgae have a different δ¹³C pattern (“fingerprint”) in AAs due to differing biosynthesis of AAs (Larsen et al., 2009; Larsen et al., 2013). It is also reported that δ¹³C patterns are largely unaffected by environmental conditions. For example, δ¹³C_AA patterns of the marine diatom Thalassiosira weissflogii did not respond to changing environmental conditions such as light, salinity, temperature and pH, despite substantial changes in bulk δ¹³C values (Larsen et al., 2015). A similar isotope pattern was reported for seagrass Posidonia oceanica and the giant kelp Macrocystis pyrifera, which showed consistent δ¹³C_AA patterns despite varying season and growth conditions (Larsen et al., 2013). The δ¹³C_AA patterns in producers, especially those of essential amino acids (δ¹³C_EAA), can be reflected in consumer tissues with little isotope effect. This is because animals obtain EAAs from their diet and EAA fractions are thought to be directly assimilated (McMahon et al., 2010). These general expectations were reasonably met in our δ¹³C_EAA dataset. Normalised δ¹³C_EAA patterns of our producer samples including mangrove leaves (yellow leaves of A. marina) and MPB did not differ between the two sites despite differing environmental conditions and substantial differences in bulk δ¹³C values (Table 4, Figs. 6 and S1). Furthermore, the consumer δ¹³C_EAA patterns also did not differ between the two sites. (Figs. 6 and S1). These findings did not support changes to feeding dependency following mangrove loss but suggested that the overall differences in the consumer bulk δ¹³C values were most likely driven by differences in the resource organic matter δ¹³C values, e.g. changes to MPB δ¹³C values that were likely associated with lower mangrove C fixation or respiratory inputs following the mangrove mortality. Furthermore, such findings indicate that the reported substantial change to the mangrove benthic faunal assemblage following the mangrove loss (Harada et al., 2019) was probably driven more by modification of physical habitat structure than changes in the use of food resources.

Figure 7 A conceptual diagram showing the ecological and biogeochemical legacy of the mangrove forest dieback in the Gulf of Carpentaria and four predicted recovery scenarios of the mangrove ecosystem with isotopic trajectories (δ represents the isotope values of animals, plants and sediment).

5 Conclusions

Reporting rare and extreme biological events can be complicated because in many cases they may occur suddenly; therefore drawing comparisons between pre- and post-event conditions remains a challenge. Our field investigations using traditional ecological techniques combined with SIA measured the initial dieback and also early recovery of an impacted mangrove ecosystem and compared an adjacent unimpacted reference system. Although the unimpacted and impacted forests were not replicated in this study, the difference between the two sites was clear. Mangrove seedling and sapling populations that increased during the period from 2016 to 2018 (8 to 32 months after the mortality event) in the impacted site suggest recovery of the mangrove vegetation. This also suggests that the environmental conditions at the impacted site are still conducive to the re-establishment of mangroves, allowing recruitment of seedlings and development of regrowth. However, mangrove leaves collected in the impacted site in 2018 showed relatively higher δ¹³C values (−25.8 ± 1.0 ‰) that are probably associated with continued water stress. Invertebrates from the impacted site representing the feeding types of grazing, leaf feeding and algae feeding were more enriched in ¹³C, ¹⁵N and ³⁴S relative to those from the unimpacted site. For example, they were more ¹³C enriched at the impacted site, by 1.7 ‰–4.1 ‰, and the difference did not change over the study period. Overall, our stable CNS isotope data supported the hypothesis that changes to biogeochemical processes occur following the mangrove mortality. These changes include lower mangrove C fixation and respiration, lower N fixation, and lower sulfate reduction. However, our isotope data did not support the second hypothesis that the isotopic compositions change over time with recovery of mangrove vegetation. Recovery of biogeochemical processes was not evident even after 2 years, suggesting an ongoing impact of the mortality event. An exception was that N cycle recovery may be occurring faster.

Considering that the environmental conditions at the site play an important role in facilitating recolonisation of mangroves, we conceptualise the recovery of the mangrove forest under four different scenarios to give insight into the ecological and biogeochemical consequences of changing forest conditions (Fig. 7): (1) the forest recovers with mangroves being able to recolonise at the site without future perturbations; (2) the forest recovers with future perturbations such as climatic events (for example, mangrove recolonisation is driven by events such as ENSO cycles); (3) the forest does not recover and is transformed into intertidal mudflats; and (4) the forest recovers partially at the site and in a reduced size and/or is recolonised by other plants such as saltmarshes (e.g. mangroves only recolonise in the lower intertidal zone). Each of these scenarios will have a distinct isotopic trajectory for C, N and S. While these are likely scenarios and there might be other possibilities, comparing the impacted forest and an adjacent unimpacted reference forest can help us quantify the recovery. Continued monitoring of the post-dieback forest would be required to predict the long-term trajectory of ecosystem recovery and how ongoing climate change and extreme climatic events affect the recovery of mangroves in the impacted region. In such a long-term investigation, SIA is a powerful tool, capable of tracking changes in biogeochemical processes over time. As such, it is of great assistance in ecosystem analyses and detecting the underlying biological mechanisms that drive changes and recovery.