Causal impact of fire on a globally rare wetland plant: a 40-year Bayesian time series analysis

Canby’s dropwort (Oxypolis canbyi (J.M. Coult. & Rose) Fernald) was listed as federally endangered in 1986, yet the species has continued to decline and is no longer found in 11 counties throughout its former range. The seasonal wetlands in which this forb occurs are disappearing from the landscape, often closing in and transitioning to wet forest or are drained and converted to agriculture. We document the effects of reintroducing fire to the only population of O. canbyi in Maryland and examine the resulting population increase using Bayesian interrupted time series analysis with a counterfactual. After cutting woody vegetation, 3.74 times more stems of dropwort per year were produced over the baseline of no intervention. Subsequently, after fire was reintroduced, 10.80 times more dropwort stems per year were produced in comparison to the time period after cutting woody vegetation. The counterfactual prediction showed that in the absence of intervention with fire, dropwort stem production would likely have declined. Cutting woody vegetation set the stage for the growth of fine fuels by increasing sun exposure on wetland grasses. It was only after fire was reintroduced that the population of plants expanded significantly in extent and flower production. A process for model selection with a directed acyclic graph followed by Bayesian interrupted time series analysis and a counterfactual was useful for causal inference. Application of fire is an important step in the recovery of the federally endangered dropwort.


Introduction
Fire regime changes over the course of at least 120 million years have shaped the evolution of plants , and fire is one of the most important factors for maintaining terrestrial biodiversity . Restoring open wetlands and their associated species can be a delicate balance between active land management and protecting sensitive habitats. Knapp et al. (2021) found that 65 plant taxa in 31 US states, the District of Columbia, and Ontario have gone extinct since European colonization due to direct anthropogenic disturbances (i.e., habitat alteration or destruction). Plant extinction has been investigated using population viability analysis and specifically in relation to fire intervals and seasonality (Quintana-Ascencio et al. 2003;Evans et al. 2010), but in some cases, these types of studies have been hampered by temporally limited datasets (Menges 2000). Reversing an extinction trend requires that we understand which activities or interventions are most effective, while addressing limitations common in ecological datasets. Our study utilizes a method to address limitations of observational data by examining two interventions, woody plant removal and application of prescribed fire on a single rare plant population in a seasonal wetland.
Canby's dropwort (synonym Canby's cowbane), Oxypolis canbyi (J.M. Coult. & Rose) Fernald, hereafter referred to as "dropwort, " is a rhizomatous perennial herb found in shallow wetlands of the mid-Atlantic Coastal Plain. The range of this federally listed species is sometimes given as Delaware to Georgia (Tucker et al. 1983), but the Delaware population is now considered to be extirpated (USFWS 2022). Hamrick et al. (2019) found the Maryland population to be genetically isolated and disjunct from the rest of the populations found further south (North Carolina is the next closest location). This plant grows in Coastal Plain habitats, including pond cypress savannas, wet meadows, and depressional wetlands. Although these habitats occur throughout its range, most do not currently support the open savannah conditions necessary for the dropwort (USFWS 2015).
Delmarva bays are seasonally flooded depression wetlands found along an Atlantic coastal peninsula shared by Delaware, Maryland, and Virginia (Stolt and Rabenhorst 1987;Phillips and Shedlock 1993;McAvoy and Clancy 1994;McAvoy and Bowman 2002). These wetlands, characterized by long periods of inundation and open canopy, have declined throughout their historical range due to drainage, conversion to farmland or pine plantations, and fire suppression (NatureServe 2022). Habitat loss and wetland degradation combined with a lack of natural or prescribed fire have resulted in a continued decline of dropwort populations since its listing in 1986 (USFWS 2022). Evidence of fire in this wetland type was found in dated sediment cores from a nearby Delmarva bay, which found a general decline in grass pollen, an increase in deciduous tree pollen, and a decrease in microscopic charcoal since European settlement (unpublished report (Sneddon 1991). In addition, an index of fire based on soil charcoal has been related to vegetation type on the Delmarva peninsula (Kirwan and Shugart 2000) and indigenous fire knowledge has been passed down to the present day in this region (J. Kirwan, personal communication, Maryland, USA, 2022). When indigenous fire stewardship of beargrass (Xerophyllum tenax [Pursh] Nutt. Melanthiaceae) was simulated, the population growth rate increased after fire and traditional harvest (Hart-Fredeluces et al. 2021). Experiments have also shown increased growth of a plant species (Anthoxanthum nitens (Weber) Y. Schouten & Veldkamp) after harvests due to increased tillering (Reid 2005) and increases in germination using water infused with smoke (Shebitz and James 2010). No similar experiments have been published about this dropwort species, as the rarity of occurrences and endangered species status makes replicated experiments impractical. In a field guide to rare plants in Georgia, Chafin (2007) recommended the application of prescribed fire at two to three 3-year intervals for the management of dropwort populations. The 1990 USFWS recovery plan only mentions fire twice, whereas the most recent 5-year review (2022) includes a recommendation to remove "shrub/tree encroachment with prescribed fire, canopy thinning, or other techniques. " (USFWS 1990(USFWS , 2022. To study methods for increasing the viability of this dropwort population, we asked the following questions: Would manual removal of encroaching woody vegetation be sufficient to perpetuate the dropwort population, or would controlled burning be more effective at increasing the population size? To address these questions, we followed a causal inference process outlined by Laubach et al. (2021) with Bayesian time series modeling. Bayesian time series analysis and the counterfactual approach are used more often in diverse fields from medical interventions to economics and marketing, but have been used recently to analyze wildlife census data (Bucciarelli et al. 2020;Jellesmark et al. 2021) and the impact of reservoir drawdown on water clarity (Nakanishi et al. 2022). The methods we present here to interpret the results of monitoring data from a single rare plant population include implementing interventions where randomized trials are not possible, followed by modeling and prediction of a counterfactual (Brodersen et al. 2015;Natesan Batley et al. 2021).

Study site
The dropwort population is located in a single Delmarva bay on the Eastern Shore of Maryland near the Delaware state line at a preserve managed by The Nature Conservancy (TNC) (Fig. 1). This preserve was acquired in 1983, although monitoring of the dropwort began in 1982, with 36 individual stems (Boone et al. 1984). The location of Soils at this location consist of Corsica mucky loam, with 0-2% slope (Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture. Web Soil Survey. Available online at http:// webso ilsur vey. nrcs. usda. gov/ accessed January 19, 2022). This is consistent with other dropwort sites, which tend to be sandy loam or loam soil with a medium to high organic content, high water table, and deep, poorly drained and acidic soils underlain by a clay layer (NatureServe 2022).

Propagation and girdling
In 1989, an attempt was made to artificially increase the dropwort population in Maryland by removing two of four remaining plants and propagating them at the North Carolina Botanical Garden. Sixty cultivated plants were subsequently re-introduced to the site in 1992 in 13 locations in the wetland and monitored annually. Transplants declined every year until no transplants were observed in 1999. Thus, propagation was not found to be a viable tool for dropwort protection.
When the preserve was acquired in 1983, woody vegetation encroachment by red maple (Acer rubrum) and sweet gum (Liquidambar styraciflua) was substantial. The trees were shading out herbaceous vegetation, with portions of the wetland completely under a tree canopy. Additionally, there was concern about increased evapotranspiration as the trees increased in size, lowering the water levels in the wetland. In 2002, the majority of trees within the wetland were removed or girdled. Cut stumps and girdled trees were treated with 50% glyphosate (cutting phase). Re-sprouting seedlings were manually removed in subsequent years.

Intervention with prescribed fire
It was determined from multiple sources (Kirkman 1995;Chafin 2007) that a controlled burn might be the most effective means of maintaining the open character of this wetland. Additionally, the native grass Panicum hemitomon had spread substantially following woody vegetation removal (Fig. 1a). This grass was encroaching on the dropwort area, creating a thick layer of thatch and possibly preventing dropwort seeds from landing on a suitable substrate for germination. Planning for a burn was problematic, however, as the sensitive nature of the wetland made the creation of standard fire breaks undesirable. Thus, a burn prescription was prepared that allowed for fire to burn off the vegetation in the open wetland but not spread into the surrounding shaded woods, utilizing the fuel type transition as a natural firebreak.
Prescribed fires were successfully implemented in October 2015 and 2017 after the dropwort had senesced.
The burn prescription required that the relative humidity and fuel moisture be high enough that fire would not carry into the woods, thus only burning the herbaceous vegetation and other fine fuels within the pond where the dropwort was present. On the day of each burn, a test fire was ignited in the woods and observed. If it self-extinguished, it was determined within the prescription and the burn proceeded. Each burn was conducted at midday and lasted approximately 1 h (burn #1: 23C, winds S 8MPH, RH 52%, burn #2: 28C, winds N 6MPH, RH 58%). Immediately post-burn, the dropwort site was inspected; on each occasion, above-ground vegetation was consumed by fire up to the base of the dropwort stems, where soil was still damp (Fig. 1b, c).

Monitoring
Monitoring was conducted in mid-August when the dropwort was in bloom (Fig. 1a), making identification of individual stems easier. This study analyzes the results of plant counts conducted from 1982 to 2021. Initial monitoring consisted of counting each individual dropwort stem, measuring its height, and recording whether it was in bloom, fruit, or vegetative. However, due to the high number of stems after 2016, it was no longer possible to count nonflowering stems without damaging the population; thus, only flowering stems were counted from 2017 to 2021.

Precipitation and temperature data
Data from nearby weather stations was obtained from the National Climate Data Center Archive (NCDC 2022). Monthly summaries of total precipitation were divided into fall, winter, spring, and summer periods and analyzed for correlations with stem count data. The Palmer Drought Severity Index (PDSI) and Standardized Precipitation-Evapotranspiration Index (SPEI) values were obtained from the North American Drought Atlas (Cook et al. 2010;Beguería et al. 2014).
After we determined important factors and confounders for the dropwort data, we constructed a directed acyclic graph (DAG, Fig. 2). A DAG is a visual representation of the study subject based on the structural causal model framework (Pearl 2009) that makes assumptions explicit and reduces bias for causal inference from observational data (Arif and Aaron MacNeil 2022). Two different multivariable model selection processes were used to account for confounders shown by the DAG, test for interactions and predict a counterfactual.

Model selection
An interrupted time series model with first-order autocorrelation (ITSA, (Natesan 2019)) was determined to be best suited for analyzing our dataset. The data are time series, taking place over a 40-year period of time with repeated annual measurements. The data were assumed to be autocorrelated because dropwort is a perennial; the number of plant stems in a given year has an impact on the number of stems the following year. The ITSA model type is appropriate because the interruptions in the time series are the introduction of interventions such as cutting and prescribed fire. The dependent variable was the stem count of dropwort. Bayesian estimation was used to provide more information about the parameters in the form of a posterior distribution. The posterior distribution contains all probable values of the parameter and therefore the credible interval estimates obtained are straightforward to interpret. More details about the methodology specifically for this kind of data can be found in Natesan (2019) and Natesan Batley (in press). The Bayesian rate ratio (BRR) effect size proposed in Natesan Batley et al. (2021) was computed to measure the magnitude of the effect of each intervention on stem count. Therefore, two BRRs were estimated. BRR 1 measured the ratio of the rate of stem count increase between the cutting phase and the baseline phase where no treatment was applied. BRR 2 measured the ratio of the rate of stem count increase between the fire phase and the cutting phase. Larger BRR values would indicate that there is a higher stem count increase in the phase under consideration than the previous phase. The method of calculating this ratio is given in the following equations: Let θ jt be the expected number of stems counted at time t in phase j where for the baseline phase, j = 0. The number of stem counts in j at time t is given as Y jt ∼ Pois θ jt .
If λ j is the stem count rate in j, and ρ the autocorrelation, In Eqs. 1 and 2, autocorrelation is applied within a given phase only past the first time point. The number of events at the proceeding time points depends on the number of events at the previous time point. The stem count rate λ j is given as where ξ is the stem count rate when j = 0 and exp(ω) is the rate ratio effect size (Kunz et al. 2015). The rate ratio effect size is computed as, (2) θ jt = j e j , if t = 1 j e j (1 − ρ) + ρY j(t−1) , otherwise .
(3) j = ξ exp ωj , When using predictor variables, Eq. 1 was expanded as: In Eq. 5, γ k refers to the kth effect with a total of K oneway and two-way interaction effects.
The priors that were used are specified below: Four models were fitted to the time series using a Poisson distribution. Poisson distribution was used to model the dependent variable due to the nature of the count data. Model 1 is a simple ITSA model estimated without any slopes, that is, time (year) was not included in the model. The interruption in the ITSA model is change in phase. Model 2 is an ITSA model with time included as a predictor. Each phase was allowed to have its own intercept and slope. In model 3, relevant predictor variables were added to model 2. These variables included October precipitation, October temperature, March temperature, the minimum temperature in March, previous year's PDSI, the interaction between all these variables and phase indicator (i.e., cutting phase, fire phase), and the interaction between October temperature and October precipitation. For model 4, variables that were not statistically significant were removed from model 3 to form a more parsimonious model. Deviance information criterion (DIC) was used to retain the models. Similar to the Akaike information criterion (AIC), lower values of DIC are used to select models which have better fit and fewer effective parameters (Spiegelhalter et al. 2002;Spiegelhalter et al. 2014).
A plant condition index was calculated using the first year of the counts as the baseline by dividing each subsequent year by the value for the first year (1982 = 1 in the index). The Bayesian rate ratio from the previously described ITSA model was used to pick the most important intervention phase (reintroduction of fire). Using the R statistical software (version 4.2.1, R core team 2022), four additional models were fitted, using the BSTS package (Scott and Varian 2013). Model 5 is a state space that includes auto-regression (AR), whereas model 6 includes AR and a student linear local trend. The state space for model 7 included AR and dynamic regression on covariates (October precipitation, October temperature, SPEI, and minimum March temperature). Model 8 state space included AR and covariates as model terms (October precipitation, October temperature, SPEI, and minimum March temperature).
The BSTS package was used to compare one ahead errors through each time step for all four models. The model with the lowest overall error (model 8) was then used in the CausalImpact package to predict a counterfactual for comparison to the actual data after the intervention (Brodersen et al. 2015). The DIC for ITSA models 1-4 respectively were 4566.9, 3099.18, 2900.16, and 2764.871. In model 4, the rate ratio between phases 1 and 2 was 3.74 [2.5, 5.31] and phases 2 and 3 was 10. 80 [8.58, 13.27], indicating that there was a statistically significant difference in the stem counts between the phases. More importantly, fire was the intervention with the highest rate ratio. Figure 3 shows the posterior distribution of the effect sizes, that is all possible values the effect size could take along with their frequency. This helps us understand not just the effect size estimate, but also how strong the evidence is. For instance, a narrow distribution would indicate that more values are surrounding the effect size estimate, indicating we have stronger evidence. In this figure, we can see that the effect size for cutting is approximately 3.74 and for fire it is 10.80.

Post
The slopes for the trend of phases 1 (no intervention), 2 (cutting), and 3 (fire) were 2.57 [1.05, 4.07], −2.96 [−6.29, 0.4], and 15.71 [10.06, 21.33], respectively. This indicates that there was some increase with time during phase 1 but statistically no increase with time for the cutting phase (because the credible interval contains zero), and a much more prolific increase in the fire phase (Fig. 4, Table 1). October precipitation (b = −6.32 [−7.79, −4.89]), previous year's PDSI (b = 5.31 [0.82,9.96 Although higher precipitation in October was associated with lower stem count, higher precipitation in October was associated with higher stem count when fire or cutting was used. Although March temperature by itself did not have a statistically significant effect on stem count, lower March temperatures were associated with higher stem count in the cutting phase while higher March temperatures meant higher stem count in the fire phase. The rest of the main and interaction effects' credible intervals contained zero and, therefore, were deemed statistically not significant (Table 1).
The counterfactual prediction for the dropwort without a fire intervention versus observed stem production using the BSTS and CausalImpact packages in R is shown in Fig. 4. During the post-fire intervention period, the dropwort index (1982=1) had an average value of 40.46. By contrast, in the absence of an intervention with fire, we would have expected an average index value of 2.46 [95% CI −0.11, 5.14]. Subtracting this counterfactual prediction from the observed response yields an estimate of the causal effect the fire intervention had on the dropwort index, 38.01 [95% CI 35.32,40.57]. Tyndall (2000) suggested that Delmarva bay wetlands are adapted to fire. However, applying fire in this wetland was delayed because of the reluctance to potentially negatively impact this sensitive site and the only population of the federally endangered dropwort in Maryland. It was not until we established a burn prescription that allowed for a minimally invasive burn that we were comfortable proceeding with this intervention.

Discussion
By examining aerial photography, Tyndall (1989, unpublished report to TNC) found that without management, herbaceous cover decreased by 83% from 1938 to 1980 in Delmarva bays, with a 58% drop from 1952 to 1980 at our study site. Extinction risk from the loss of discrete small seasonal wetlands can be disproportionately higher than from a similar loss of extent in continuous habitats (Deane et al. 2017). Conceptual models have been developed from data collected in Carolina Bays that include both hydrology and fire frequency as important factors for maintaining herbaceous diversity, similar to our findings (Kirkman 1995;Kirkman et al. 2000). The results of overstory tree removal initiated in 2002 were positive, with a small step increase in dropwort numbers but no significant growth trend over the subsequent 14 years (Table 1, slope of phase 2). Dropwort was potentially being suppressed by P. hemitomon, which was responding aggressively to the increased light availability, accumulating a thick layer of thatch. Panicum hemitomon responds with increased growth after fire, but is delayed in emergence, and spring growth is reduced after a fire followed by inundation (Kirkman andSharitz 1993, 1994). In this region, geographically isolated wetlands' hydrology is influenced by a wet season that is energy limited from June to November and a waterlimited dry season from December to May ) (see Fig. 1a-e for the dry season and Fig. 1f for the end of the wet season). During a 2013 survey, it was observed that the P. hemitomon thatch appeared to be preventing dropwort seed germination outside the small area in which it occurred (W. Knapp, personal communication, MD-DNR Heritage Service, Maryland). In a survey of dropwort populations in Georgia, Bowling (Bowling 1986, The Status of Canby's Dropwort (Oxypolis canbyi) in Georgia. Unpublished report submitted to the USFWS by the TNC Eastern Heritage Taskforce) observed that native Carex and Panicum species often exclude dropwort by growing too densely. Despite the observed small stepwise increase, our counterfactual model showed that woody vegetation removal alone would likely have resulted in a slight decline of dropwort stems over time.
The controlled burn did result in the removal of thatch, creating large patches of exposed soil immediately postburn (Fig. 1c-e). Research in a greenhouse experiment found that the removal of previous years' growth results in increased light availability which favors brackish marsh plants able to take advantage of the increase (Bickford et al. 2012). Increased light availability after the prescribed fires could help explain our post-fire observations of 10.80 times more dropwort stems per year, in comparison to the time period after only cutting woody vegetation (Fig. 3).
The counterfactual prediction showed that in the absence of the intervention using prescribed fire, dropwort stem production would likely have declined given the trends in weather variables following 2016 (Fig. 4). Future research might investigate other fire effects on the wetland, such as nutrient cycling, pyrogenic carbon, and light availability at ground level, to help establish more specific causal factors from fire effects leading to increased dropwort expansion. A counterfactual approach to measure impact will continue to be helpful in guiding conservation interventions for this species and could be applied more broadly to other species of concern or other rare habitat types (Grace et al. 2021).

Conclusions
Controlled burning was the most effective tool for increasing the Maryland dropwort population. While woody vegetation removal did result in an increase in population size, it was not until the wetland was burned that the population increased significantly.
Additionally, after fire, the dropwort population expanded beyond the extent of the very small patch in which it had historically occurred, and the percentage of stems flowering increased.
This study highlights the challenges land managers are faced with when deciding on the most appropriate restoration tool. In a fire-suppressed environment, fire surrogates such as the manual removal of woody vegetation can result in positive effects. However, unforeseen consequences of restoration activities (e.g., increased competition) and potential negative impacts (e.g., fireline creation) must all be taken into consideration when selecting management practices that provide the greatest restoration outcomes.