The role of understory phenology and productivity in the carbon dynamics of longleaf pine savannas

Savanna ecosystems contribute ~30% of global net primary production (NPP), but they vary substantially in composition and function, specifically in the understory, which can result in complex responses to environmental fluctuations. We tested how understory phenology and its contribution to ecosystem productivity within a longleaf pine ecosystem varied at two ends of a soil moisture gradient (mesic and xeric). We used the Normalized Difference Vegetation Index (NDVI) of the understory and ecosystem productivity estimates from eddy covariance systems to understand how variation in the understory affected overall ecosystem recovery from disturbances (drought and fire). We found that the mesic site recovered more rapidly from the disturbance of fire, compared to the xeric site, indicated by a faster increase in NDVI. During drought, understory NDVI at the xeric site decreased less compared to the mesic site, suggesting adaptation to lower soil moisture conditions. Our results also show large variation within savanna ecosystems in the contribution of the understory to ecosystem productivity and recovery, highlighting the critical need to further subcategorize global savanna ecosystems by their structural features, to accurately predict their contribution to global estimates of NPP.


INTRODUCTION
Understory vegetation is thought to play a key role in carbon dynamics of ecosystems across the globe. Uncertainties exist, however, regarding the total carbon found in this vegetation strata (Johnson et al. 2017, Refsland andFraterrigo 2017), and these uncertainties are expected to increase as climatic extremes become more common and severe (IPCC 2014). These shifts in climate, particularly drought, will have a direct effect on ecosystem structure (both understory and overstory), function (phenology and physiology), and ultimately their management strategies (Rich et al. 2008, Becknell et al. 2015. Quantifying phenological changes and understory carbon dynamics is particularly difficult using coarse remote sensing techniques due to the canopy variation in forests and/or uneven sensor obstruction (Picotte and Robertson 2011). Furthermore, critical variation of the understory occurs at fine (1-10 m 2 ) spatial scales (Hiers et al. 2009, Starr et al. 2015. This is especially the case for savanna and woodland ecosystems where canopies can be relatively open and less uniform, and therefore harder to classify using coarse spatial scales of satellite-derived remote sensing products (Patenaude et al. 2004, Simard et al. 2011, Smith et al. 2014. Understanding savanna structure and function is critically important, as savannas cover~15% of the global land surface (Ma et al. 2013) and con-tribute~30% to global net primary production (NPP, Field et al. 1998, Fei et al. 2017, thus representing a significant portion of the global carbon cycle (Bombelli et al. 2009). Nevertheless, savanna ecosystems vary across the globe based on their underlying geology and regional climates (Bond et al. 2003, Vaughn et al. 2015. Differences among and within savannas often result in complex responses to environmental fluctuations (Boke-Ol en et al. 2016, Starr et al. 2016. For example, savanna ecosystems were shown to exhibit different responses to drought, even within the same region (Starr et al. 2016). This finding has been linked to variations in soil waterholding capacity, which in turn has led to variations in plant functional types among sites (Vaughn et al. 2015). Differences in plant functional types can drive changes in phenological patterns, as well as differential response to disturbances (Rich et al. 2008, Bernhardt-R€ omermann et al. 2011, Wu et al. 2016. Thus, inclusion of sitespecific phenology is extremely important in making accurate predictions regarding physiological patterns and ecosystem carbon dynamics. These patterns, however, are poorly represented in global dynamic vegetation models and are a source of systematic errors in regional to global scale model predictions (Boke-Ol en et al. 2016).
Complicating our understanding of ecosystem carbon and energy dynamics is the fact that disturbance events result in differing responses between the understory and the overstory (Rich et al. 2008). The understory in forests and woodlands can contribute roughly 30-50% to the annual carbon fluxes of an ecosystem (Powell et al. 2008, Bracho et al. 2012). However, this response varies across the landscape and may change with increasing frequencies of the disturbance (Mitchell et al. 2014). A study on North American savannas revealed that the understory was more capable of rapidly utilizing newly available resources following a prolonged drought compared to the overstory (Rich et al. 2008). This rapid shift to resource availability offset the whole ecosystem decline in function and accelerated recovery from the disturbance (Rich et al. 2008). In contrast, Ma et al. (2013) showed that increasing temporal variability in phenology coincided with decreasing tree-grass ratios, which suggest that higher woody fractions in the understory promote higher stability and resilience to future climate change. Collectively, these studies highlight the influence of site-specific abiotic and biotic factors that drive savanna overstory and understory phenological patterns and ultimately the ecosystem's carbon dynamics.
Savanna understory vegetation also plays a key role in enabling the efficient spread of surface fires through these ecosystems, which maintains the ecosystem's structure and function (Barlow et al. 2015). With a bare or patchy understory, fire would not effectively spread across the floor and therefore enable the expansion of woody shrubs , leading to a shift in the ecosystem structure (Peterson and Reich 2007), and hence a shift in ecosystem carbon dynamics. Changes in structure can also precipitate declines in biodiversity (Walker and Peet 1984) which further alter the carbon cycling in savanna ecosystems (Mitchell et al. 1999, Kirkman et al. 2016a). However, the degree to which fires interact with understory regrowth, especially during drought, has yet to be fully understood. How fire alters the overall carbon fluxes of the ecosystem, as a function of understory regrowth, varies substantially across different ecosystems. A study conducted in a longleaf pine ecosystem found that physiological activity and hence productivity returned to pre-fire levels only 30-60 d post-fire (Starr et al. 2015). A separate study in a mixed longleaf slash pine forest in north central Florida (USA) suggested that 70% of understory N and C pools are consumed by managed burns, but recover within three years post-fire (Lavoie et al. 2010). These studies, however, were short in nature-following a single fire cycle, showing the need for longer term observations to cover greater environmental conditions and multiple fire cycles.
Longleaf pine savannas of the southeastern United States are ideal for these studies since their canopies are open and dominated by longleaf pine, Pinus palustris (Mill.), while their understories host a remarkably high quantity of species when compared to other temperate ecosystems Peet 1984, Kirkman et al. 2004). Up to 50 species per m 2 can be found within the understory, mostly on poorly drained soil Peet 1984, Kirkman et al. 2001), whereas drier and sandier soils often have lower biodiversity in the understory (Goebel et al. 2001. The relative species diversity in longleaf pine savannas is often attributed to differences in nutrient availabilities in the soil, as more mesic sites usually have higher proportions of clay in the soil, which has been shown to drastically lower nutrient availability (Metz 1952, Christensen 1977, Pellegrini 2016. Quantifying understory phenology, as a function of differences in understory composition, and its contribution to regional and global carbon budgets can be tedious and time-consuming and may even be problematic (Johnson et al. 2017). However, the Normalized Difference Vegetation Index (NDVI; Boke-Ol en et al. 2016) has been demonstrated to be advantageous in determining the magnitude of ecosystem carbon and energy fluxes (Potter 1999), as well as in quantifying understory phenology (Glenn et al. 2008). For example, NDVI has been used to extract phenophases (greenup, maturity, senescence, and dormancy; Zhang et al. 2003) and predict leaf phenological patterns (Xu et al. 2014) for temperate deciduous trees. NDVI was also successfully used to quantify phenology patterns at FLUXNET sites; however, the accuracy depended on the plant functional types present at each site (Balzarolo et al. 2016). Other studies have used NDVI estimates to quantify food availability and nutritional quality for avian herbivores, showing that using handheld instruments to measure NDVI had a higher accuracy compared to satellite-derived measures (Hogrefe et al. 2017).
The spatial resolution of NDVI measures is critical for interpretation of ecological relevance. Studies have found that utilizing temporal and/or spatial resolutions that are too coarse lowers prediction accuracy when studying interannual phenological patterns or biomass productivity (Gatis et al. 2017). The course nature of satellite-derived NDVI measures in evergreen forests has been shown to be a better measure of ecosystem structure and function, rather than gross ecosystem productivity (GEP; Restrepo-Coupe et al. 2016). These studies highlighted the usefulness of NDVI measures and its limitations with respect to coarse spatial scales. This underlines the need for increased use of ground-based remote sensing techniques, especially when site-specific characteristics play an important role in the productivity of the ecosystem of interest (Hogrefe et al. 2017).
Using sub-canopy remote sensed products, the overall objective of this study was to determine the contribution of understory phenology to ecosystem CO 2 fluxes at the ends of a soil moisture gradient with respect to prescribed fire and seasonal drought. We used monthly understory NDVI measurements to quantify productivity at two longleaf pine sites representing the ends of this gradient (mesic and xeric). We hypothesized that (1) the xeric site would recover more rapidly from drought as it evolved under lower soil moisture conditions; (2) fire would lower understory NDVI at both sites, with the mesic site understory recovering more rapidly post-fire due to higher soil water availability; and (3) understory NDVI would be associated with increased ecosystem productivity at both sites, with the xeric site exhibiting greater influence of the understory on overall carbon fluxes due to its more open canopy.

Study site
The study was conducted from August 2011 through to March 2016 in southwest Georgia, USA (31.2201°N, 84.4792°W), within the Southeastern Coastal Plain, at the Joseph W. Jones Ecological Research Center (JWJERC; Fig. 1). The area receives on average 1310 mm of rain annually and is climatically categorized as humid subtropical (Christensen 1981, Goebel et al. 2001). The long-term average monthly air temperatures range from 9.6°C to 27.6°C in January and July, respectively.
Locations for this study were established in two sites that represent the xeric and mesic ends of an edaphic gradient in longleaf pine ecosystems at the JWJERC and are located~5.0 km from one another ( Fig. 1; Kirkman et al. 2001). The xeric site was characterized by deep sandy soils classified as typic quartzipsamments without an argillic horizon in the upper 300 cm of soil, with a water-holding capacity of approximately 18 cm of water per meter of soil, because of the high porosity and underground drainage in the soil (Goebel et al. 2001). The mesic site had a soil composition that consisted of sandy loam that rested on top of either clay or sandy clay loam soils. The mesic  site soils were not well drained and characterized as arenic paleudults. The argillic horizon of these soils was within~95 cm, and the water-holding capacity was 40 cm per m of soil in the upper 300 cm (Goebel et al. 2001. The relatively small discrepancies between the site characteristics contributed to a large distinction between the sites' species composition and vegetation structure ). Both sites were comprised of an open longleaf pine canopy in the overstory and wiregrass (Aristida beyrichana Trin.) in the understory, while the composition of secondary species varied by site . Both ecosystems were woodland systems with basal areas of 11.1 and 18.4 m 2 /ha at the xeric and mesic sites, respectively, of which~18% were oak trees at the xeric site versus 6% at the mesic site. Due to differences in soil drainage, understory species richness was significantly higher at the mesic site compared to the xeric site (Kirkman et al. 2016a). The understory of the xeric site, however, had a higher density of native hardwoods relative to the mesic site, with the scrub oak species Quercus laevis Walt. and Quercus margaretta Ashe found most prevalently. The mesic site had a lower concentration of oaks present, and instead, a higher population of common persimmon Diospyros virginiana L. found most frequently in the understory (Goebel et al. 2001). The overstory leaf area index (LAI) for the mesic site ranged from 0.65 to 1.1 m 2 /m 2 , while at the xeric site, it ranged from 0.22 to 0.65 m 2 /m 2 (Addington et al. 2006, Wright et al. 2013). Both the mesic and xeric sites have been thoroughly described in previously published works by Mitchell et al. (1999), Kirkman et al. (2001), Ford et al. (2008), and Starr et al. (2015).
Four productivity plots, 50 9 50 m in size, were established prior to this study at each site. These plots were located within the flux footprint of each site's eddy covariance tower (~500 m radius from the tower). Understory composition at both sites was estimated annually during each fall within the four plots, with the exception of 2014 and 2016, due to a change in the sampling protocol from annual measurements to 2-yr measurement cycles in 2013. The understory biomass was estimated using 0.75-m 2 clip plots. A plot frame was randomly tossed from each of seven litter trap positions, which had been previously established within the plots at the two sites, resulting in a total of 28 replicates at each site. Where the frames touched the forest floor, all live and dead vegetation smaller than 1 meter in height was clipped and brought to the laboratory for analyses. The vegetation was then classified by growth form into forbs, ferns, legumes, other grasses, and woody plants. The biomass from the plots was then dried to a constant weight.

Drought indices
The Palmer Drought Severity Index (PDSI) for the study period was obtained from the National Oceanic Atmospheric Administration (NOAA) and the National Centers for Environmental Information (NCEI), Ashville, North Carolina, USA (National Oceanic and Atmospheric Administration, Department of Commerce n.d.). The index describes severity of a drought on a scale from À5 to +5, where more negative numbers indicate more severe drought. The index is computed using air temperature, rainfall, and soil moisture measurements at climate stations around the globe (Palmer 1965, Dai et al. 2004).

Prescribed burns
The two sites were burned frequently prior to the study, averaging a two-year fire return interval since 1994 with the last burn prior to the study occurring in 2009. During our study, the sites experienced 3 burns, during springs of 2011, 2013, and 2015. Strip head fires were set every 30-50 m (depending on the local fire conditions at the time of the burn) starting downwind of the management unit and moving upwind (Hiers et al. 2009). No damage to overstory trees was recorded following the prescribed fires, and as fires typical to this surface fire regime are low intensity , flame heights are typically below 2 m, thus unable to reach overstory crowns which are on average 23 m above the forest floor (Robertson and Ostertag 2007, Varner et al. 2007, Whelan et al. 2013). While the two study sites had prescribed fires on a twoyear cycle, there was no control site that was absent of prescribed fires, due to rapid ecosystem movement toward an alternative stable state when deprived of prescribed fires (O'Brien et al. 2010). If prescribed fire is withheld for as little as 4 yr, the ecosystem structure, composition, and function change (Kirkman et al. 2013).
In addition to data collected describing understory vegetation composition, we sampled aboveground litter and understory biomass before and after the fires, using methods described in Ottmar et al. (2007). The number of clip plots sampled at the two sites varied from 10 to 20, such that a standard error of the mean of 15% was achieved (Whelan et al. 2013). In 2011, 2013, and 2015, we collected aboveground biomass within two weeks pre-and post-fire using paired 0.75-m 2 clip plots at the mesic and xeric sites (Table 1). Clip plots at both sites were located along a transect and sampled every 25 m starting at the base of each flux tower and extending within the flux footprint (~500 m) in the direction of the prevailing wind (Whelan et al. 2013). Harvested clip plot litter and biomass were dried to a constant weight and mass, and fuel consumption was calculated as the difference of pre-burn and post-burn dry weight (Whelan et al. 2013).

NDVI measurements
Understory NDVI was measured monthly using a Tetracam ADC camera attached to a tripod, which was placed consistently at 1.60 m above each subplot for each measurement. Because the xeric site was spatially more complex than the mesic site, four plots were established within the footprint at the xeric site, whereas the mesic site was sufficiently represented by only three plots. In both sites, four 1m 2 subplots were located 10 m from the center of each plot, one in each of the four cardinal directions (north, east, south, and west). Because the camera field of view could not span the entire 1 m 2 , four pictures were taken in each of the four subplot corners, which were processed electronically and averaged during the post-processing procedure. Before each subplot measurement, the camera was calibrated using a white Teflon calibration chip, to account for changes in reflective values due to weather conditions. Post-processing was performed using the Pixelwrench II software (Tetracam, Chatsworth, California, USA), which calculates understory NDVI from the red (Red) and near-infrared (NIR) reflectance values of each image as NDVI ¼ ðNIR À RedÞ=ðNIR þ RedÞ:

Ecosystem flux measurements
Eddy covariance towers, positioned at each site, measured net ecosystem exchange of CO 2 (NEE), as well as meteorological variables such as rainfall, air temperature (T air ), and photosynthetic active radiation (PAR) above the canopy (Whelan et al. 2013). In addition, soil moisture (SWC) and soil temperature, as well as ground heat flux, were measured in a soil array adjacent to the tower. NEE was measured via an open-path eddy covariance method through a simplification of the continuity equation and converted from lmolÁm À2 Ás À1 to g C/m 2 per 30 min. Missing half-hourly data were gap-filled using separate functions for NEE during daytime and nighttime. When photosynthetically active radiation (PAR) was ≥10 lmolÁm À2 Ás À1 , daytime NEE data were gap-filled using a Michaelis-Menten approach, and when PAR was <10 lmolÁm À2 Ás À1 , nighttime NEE data were gap-filled using a modification of the Lloyd and Taylor (1994) approach, both on a monthly basis. Where too few observations were available to produce stable and biologically reasonable parameter estimates, annual equations were used to gap-fill data by site. Gross ecosystem exchange (GEE) was calculated from the difference of NEE and ecosystem respiration (R eco ) as GEE = ÀNEE + R eco . NEE, GEE, and R eco were then summed up to monthly estimates. Daytime R eco was estimated using nighttime gap-filling equation to partition NEE. A detailed explanation of the equipment and flux processing method can be found in Starr et al. (2016) and Whelan et al. (2013). Finally, we converted NEE and GEE to net ecosystem productivity (NEP) and gross ecosystem productivity (GEP), where positive values

Statistical analysis
Prior to the analysis of NDVI and ecosystem fluxes, we formulated simple linear mixed models to quantify differences in SWC, T air , and PAR between the two sites, using site, year, fire-cycle time (FCT), and their interactions as independent variables. The factor FCT was included to describe the amount of time that had passed since prescribed burning. Following Whelan et al. (2013), FCT 1, 2, and 3 corresponded to the first, second, and third 28 d following fire, respectively. These categories comprised spring and the beginning of summer. FCT 4 and 5 represented the next 140 and further 225 d following FCT 3, encompassing late summer and fall. Pre-fire (pf) indicated the time before the next fire (~10 months). We used the lme function from the package nlme in R version 3.4.3 (R Core Team 2013), including a random effect to account for monthly repeated measurements and the subplot design, as well as an autocorrelations structure to account for the correlation between adjacent time points.
We then estimated linear mixed models for each of the response variables: monthly mean NDVI and monthly sums of NEP, GEP, and R eco . For the NDVI model, we included fixed effects for site, PDSI, FCT, mean monthly PAR, and mean SWC, as well as the interactions of site with PDSI and FCT, FCT and PAR, and FCT. The NEP, GEP, and R eco models included fixed effects for year, month, FCT, site, PDSI, and the interactions of site with FCT, PDSI, and NDVI, as well as the interaction between FCT and NDVI. Mean PAR and SWC were also included as fixed effects in these models. We also tested the models using T air instead of PAR, which resulted in similar trends of the response variables, as well as similar AIC (Akaike information criterion) values. We could not include both variables as their high correlation (>0.8) would have led to multicollinearity in the independent variables. We chose PAR for our final models to quantify changes in NDVI and ecosystem fluxes in response to changes in incoming energy. For each response variable, we used a modified backward selection criterion to choose the best, most parsimonious model. Non-significant (P > 0.05) variables were eliminated, as long as it resulted in a lower (better) AIC value. Where interactions were significant, the underlying simple effects were retained in the model regardless of significance. We examined the nature of significant differences of all interactions using the lsmeans function in R (Lenth 2016), whereby marginal means are computed holding all other independent variables at their average values.

RESULTS
Annual rainfall was~800 mm during 2011 and 2012 ( Fig. 2A), lowering the PDSI ( Fig. 2A) to <À4. Annual rainfall sums increased tõ 1500 mm during 2013, which lead to an increase in PDSI to >0, followed by another decrease during 2015 ( Fig. 2A). Soil moisture was significantly lower at the xeric site compared to the mesic site for all years and levels of FCT. Additionally, SWC was significantly lower at both sites during the drought (Fig. 2B), while SWC returned to pre-drought levels at the end of 2012 as rainfall recovered ( Fig. 2A). Mean monthly PAR decreased during the wetter years, but this difference was not significant (Fig. 2C). Average PAR was significantly higher at the xeric site (by~30-60 lmolÁm À2 Ás À1 ) compared to the mesic site for all years except 2015, when PAR was significantly higher at the mesic site (by 15 lmolÁm À2 Ás À1 ). Furthermore, PAR was significantly higher at the xeric site compared to the mesic site for all levels of FCT, except for pre-fire. Values of T air were higher (>30°C) during 2011 and then decreased below 30°C during summer of the following years (Fig. 2D). Understory NDVI during the growing season in 2011 was lower at both sites with gradual increases in subsequent years (Fig. 3A) and further increases in the summer months of 2014 and 2015. During winter months, NEP and GEP were lower at the xeric site compared to the mesic site (Fig. 2B, C). R eco decreased slightly at both sites during 2015, but GEP was only slightly lower during the winter months of that year (Fig. 3D).

Understory NDVI phenology under drought and fire
Forbs, grasses, and wiregrass comprised a greater proportion of the understory at the mesic site, specifically outside of drought years, ❖ www.esajournals.org whereas woody understory species were more abundant at the xeric site, especially during the drought (Fig. 4).
Model results indicated that NDVI varied significantly by site, but its effect depended on SWC, PAR, PDSI, and time since fire (FCT). Pre-fire (pf) understory NDVI was slightly higher at the mesic site (~0.48), compared to the xeric site (~0.43, Table 2, Fig. 5). Following the prescribed fire, NDVI at the xeric site decreased during FCT 1 (~0.26). However, no significant change was observed from FCT 1 through FCT 3 (28-54 d post-fire, Fig. 5A). NDVI at the mesic site also decreased to~0.41-0.43 during FCT 1 and FCT 2, followed by an increase in FCT 3 to~0.47 and a further increase during FCT 4 to 0.66, which was significantly different from NDVI during FCT 2 (Fig. 5). During FCT 5, NDVI at the mesic site significantly decreased to 0.53. NDVI at the xeric site significantly increased during FCT 5 (~0.66), followed by a significant decrease in FCT 5 to~0.42. NDVI significantly decreased with an increase in PDSI at both sites, with less of a decrease at the xeric site (Fig. 5B). In addition, the xeric site had significantly lower NDVI for all levels of PDSI. NDVI significantly decreased for the levels of FCT 5 and pre-fire but showed no significant change for FCT 1-3, when average PAR increased to 500 lmolÁm À2 Ás À1 . For the interaction of FCT 4 and PAR, NDVI significantly decreased by~0.5 when PAR increased, but exceeded 1, when PAR was below 100 lmolÁm À2 Ás À1 (Fig. 5). NDVI decreased by 0.3 with increases in SWC (from 17% to 24%), but these effects were not significantly different by site.

Understory NDVI effects on NEP
Net ecosystem productivity was significantly different by site and was driven by NDVI, FCT, PAR, time since fire, and PDSI (Table 3). The mesic site NEP significantly decreased with increases in NDVI compared to the xeric site    Table 3). Monthly NEP was significantly lower at the xeric compared to the mesic site pre-fire but increased (greater carbon uptake) following the fire (Fig. 6B). NEP at the mesic site significantly increased following the fires during FCT 1 and 2 (up to 40 g C/m 2 per month). Following FCT 2, NEP decreased at the mesic site to values similar to those of the xeric site (20-10 g C/m 2 per month). Both sites were weaker carbon sinks when drought increased (Fig. 6C). NEP significantly increased when PDSI moved toward zero, with the mesic site being a greater carbon sink than the xeric site. As NDVI increased, so did NEP during FCT 2 and FCT 3, which was statistically significant. During FCT 4, the opposite pattern was observed, where higher NDVI resulted in lower NEP (Fig. 6D). Photosynthetically active radiation significantly increased the carbon sink strength at both sites, with no significant difference between the mesic and xeric sites (Fig. 7A).

Understory NDVI effects on GEP
Similar to that of NEP, there was a significant difference in GEP by site which was correlated with NDVI, FCT, and PDSI (Table 3). There was also a significant effect of NDVI with FCT and a simple effect of PAR, which was not significantly different by site (Table 3). When NDVI was below 0.5, monthly GEP at the mesic site was significantly higher (160 g C/m 2 per month) than the xeric site (145 g C/m 2 per month). When NDVI (>0.75) increased, GEP decreased slightly at the mesic site, but increased at the xeric site such that it had higher productivity than the mesic site (Fig. 8A).
Preceding the prescribed burn GEP was not significantly different at the xeric site, and this persisted through FCT 2, although NDVI was lower (Fig. 8B). During FCT 3, GEP increased from 140 to 160 g C/m 2 per month at the xeric site, which then remained at that rate for FCT 4 and FCT 5. The mesic site experienced a slight decrease in GEP during FCT 1, followed by an increase to approximately 170 g C/m 2 per month and a steep decrease in GEP during FCT 3 and FCT 4 (140 g C/m 2 per month). Then, in FCT 5, the mesic site increased its productivity to roughly 155 g C/m 2 per month. Higher NDVI significantly increased GEP  during pre-fire, and FCT 2 and 3, but decreased its magnitude during FCT 4 across both sites (Fig. 8D). Interestingly, GEP increased with decreasing PAR and PDSI at both sites (Figs. 7B and 8C).

Understory NDVI effects on R eco
Model results indicated that R eco was significantly different by site and also depended on PDSI, FCT, and NDVI. Increases in understory NDVI increased R eco at the mesic site, but only slightly changed the magnitude of R eco at the xeric site. These changes in R eco , however, were not significant across NDVI levels (Fig. 9A). R eco was significantly lower at the mesic site for all levels of NDVI, compared to the xeric site.
Fire affected R eco differently by site (Fig. 9B). R eco significantly decreased at the mesic site, following the fire during FCT 1 (122 g C/m 2 per month), slightly increased during FCT 2 (127 g C/m 2 per month) and decreased further during FCT 3 and 4 (~120 g C/m 2 per month). R eco at the xeric site decreased only slightly during FCT 1-2, followed by an increase to~145 g C/m 2 per month during FCT 3 and 4. The mesic site increased its R eco during FCT 5, back to~135 g C/m 2 per month, whereas the xeric site decreased R eco to~140 C/m 2 per month. Negative values of PDSI increased R eco at both sites. However, this was more pronounced at the xeric site (Fig. 9C). For levels of PDSI >3, both sites had monthly R eco rates of~120 g C/m 2 per month, with significantly higher rates at the xeric site throughout this study. Increases in PAR significantly decreased the magnitude of R eco , with no significant difference by site (Fig. 7C).

DISCUSSION
Our study demonstrates that understory vegetation in longleaf pine ecosystems plays a complex role in patterns of carbon uptake and release. The understory phenology and composition contribute significantly to the overall productivity of this ecosystem. More importantly, our results reveal that site characteristics play a critical role in how understories mediate ecosystem recovery from drought and prescribed fire.
We show that changes in understory biomass (higher NDVI) significantly affected GEP during recovery from fire. We found increased photosynthetic capacity between two and three months after the fires when understory biomass was high. However, this relationship shifted after four months, resulting in a decrease in photosynthetic activity when NDVI was greater (Jurik 1986, Hikosaka 2005. This altered response highlights the need to incorporate understory phenology when studying ecosystem recovery and resilience in open forest stands. It further demonstrates that the capacity of ecosystems to sequester carbon depends on changes in understory biomass over the course of the year (Chen et al. 2015).
For example, minor decreases in GEP at the mesic site occurred with increasing NDVI indicated that the understory contributed less to GEP and that the overstory was the primary driver of ecosystem productivity at the site (Yang et al. Notes: df, degrees of freedom (numerator). P-value corresponds to chi-squared statistic with df. 2017). This resulted in lower photosynthetic capacity at the site when NDVI was high, compared to the xeric site, as it is dominated by longleaf pine and grass species in the understory (~30%; Wiesner et al. 2018). Furthermore, biomass and leaf area were higher at the mesic site, compared to the xeric site, suggesting less solar radiation was reaching the understory, which in turn could decrease productivity relative to the xeric site (Starr et al. 2016). Even  though basal area, leaf area, and understory NDVI were lower at the xeric site, R eco was higher with lower NDVI. These findings were likely a result of greater maintenance respiration of the oak dominated overstory during dormancy (Amthor 1984, Wu et al. 2014. In addition, the understory at the xeric site had a proportionally higher contribution to the flux of  GEP, compared to the mesic site. This effect was a result of higher abundance of woody species in the understory and therefore higher photosynthetic activity during peak NDVI, as broadleaved shrubs and trees are known to have higher photosynthetic capacity compared to grasses and coniferous trees (Weaver and Mogensen 1919, Klein et al. 2013, Renninger et al. 2015. We show that longleaf pine savannas are resilient to the disturbance of fire and drought, as a result of their ability to rapidly mobilize nonstructural carbohydrates when photosynthate supply is low, in agreement with a previous study conducted at our sites (Aubrey and Teskey 2018). However, we found that the magnitude and timescale of recovery was site-dependent. The mesic site recovered more rapidly from the disturbance of fire, possibly as a result of a higher species diversity in the understory (Kirkman et al. 2016b), which increases the potential of certain understory species to rescue ecosystem function (Elmqvist et al. 2003, Starr et al. 2015, Wiesner et al. 2018. Understory NDVI and NEP at the mesic site increased substantially more during the first three months following the fire, indicating rapid plant regrowth in the understory, potentially as a result of more readily available nutrients following the fire (Lavoie et al. 2010). This is supported by a previous study conducted at these sites, which found increased soil respiration primarily at the mesic site for two months after prescribed fire. Our findings and the findings from Wiesner et al. (2018) suggest that more rapid understory plant biomass regrowth at the mesic site increased respiration rates, as a result of higher maintenance and growth respiration. The grass dominated understory at the site also likely accelerated its recovery time, as grasses make use of stored and available resources more rapidly compared to woody species (Brewer 2011). Furthermore, higher soil moisture conditions at the site enabled faster resource acquisition (White et al. 2016), which enhanced plant biomass regrowth in contrast to the xeric site.
On the contrary, lower understory NDVI at the xeric site following the fire was likely a function of lower fuel consumption compared to the mesic site. Due to higher heterogeneity of the understory vegetation on xeric sites, fire did not spread as uniformly (Duncan et al. 2015). This would also contribute to smaller magnitude in biomass regrowth following the fires (Pinno and Errington 2016), as reflected in the ecosystem fluxes compared to the mesic site. Even though average N mineralization levels in the upper soils were found to be higher at xeric sites (Boring et al. 2004), the heterogeneity in fire effects may have resulted in lower nutrient availability within the upper soil layers post-fire, compared to the mesic site. Furthermore, similar N mineralization rates at mesic and xeric sites below the upper soil surface (<15 cm), as well as higher SWC and higher root biomass at the mesic site (Hendricks et al. 2006), would enable greater nutrient availability and uptake (Lavoie et al. 2010). Moreover, sandy soils at the xeric site can promote leaching of nutrients (Nguyen and Marschner 2013), thereby decreasing nutrient availability and thus plant regrowth.
However, the xeric site experienced a substantial increase in NDVI, GEP, and R eco 4-8 months following the prescribed burns. High R eco and GEP at the xeric site were indicative of the seasonal increase in greenness and therefore metabolic activity of deciduous oaks (Jurik 1986, Hikosaka 2005, Powell et al. 2005, thereby increasing growth respiration, as more recent photosynthates were invested into bud and leaf expansion (Keel and Sch€ adel 2010, Alla et al. 2013, Kuster et al. 2014, Herrmann et al. 2015. In support of this notion, we found a decrease in NEP for increasing NDVI during that time, suggesting that more plant biomass in the understory increased respiration. In contrast, three months following the burn, R eco was lower at the mesic site, but GEP continued to decrease for FCT 4, which also resulted in lower carbon sink capacity at the site. A similar response for the mesic and xeric sites was found in Whelan et al. (2013), who showed that the base respiration rate decreased at the mesic site but increased at the xeric site following the fires. Another explanation for lower metabolic activity at the mesic site during FCT 4 and 5 was that average T air was higher (>20°C) compared to T air for FCT 1-3 (<20°C, evident from Fig. 2D), which would cause stomata to close to conserve water (Gonzalez-Benecke et al. 2011, Samuelson et al. 2012. For example, longleaf pine trees were shown to decrease their growth rates in response to warmer summer temperatures (Foster and Brooks 2001), whereas oak species are more resistant to higher T air (i.e., anisohydric response; Roman et al. 2015).
Understory species at the xeric site were shown to be less drought resilient compared to those of its overstory. Even though the xeric site is chronically exposed to lower soil moisture conditions, we found that understory NDVI was significantly lower for all levels of PDSI compared to the mesic site, showing that understory species were not able to thrive at the level of that of the mesic site under drought conditions. However, the magnitude of NDVI change (i.e., slope of the regression) was lower at the xeric site when severity of drought increased, indicating that plant functional types at the site could buffer drought conditions more effectively (Coble et al. 2017, Granda et al. 2018, such that understory NDVI decreased by half as much as at the mesic site. In addition, overall ecosystem productivity at xeric sites was shown to recover more rapidly from drought (Starr et al. 2016), compared to more mesic sites, which was the result of drought-resistant plant species in the over-and understory (Wright et al. 2013). Oak components have been shown to increase climate resilience on xeric sites by facilitating longleaf pine regeneration . We also showed that differences in site productivity during drought were mainly governed by differences in photosynthetic activity, such that the xeric site had similar levels of R eco over the range of PDSI but could not maintain similar levels of GEP, compared to the mesic site.
In this study, the variation in phenology of the understory vegetation played an important role in the ecosystem's response to disturbance. The effects of prescribed fire on NEP and post-fire recovery were mediated not only by site type but also by interactions with variation in climate (PDSI). The ability to sustain relatively high photosynthetic capacity (and respiration) during drought at both sites has profound consequences in relation to future changes in climate, if droughts become more severe and frequent (IPCC 2014). However, prolonged drought could change the system from a carbon sink to a source, as R eco would increase over GEP, which may lead to carbon starvation or hydraulic failure, depending on the severity of the event (McDowell 2011, Zeppel et al. 2011, Klein 2015. Furthermore, severe drought could cause a shift in ecosystem structure, as a consequence of the inability to manage these systems with prescribed fire, as drier conditions limit prescription options for managers (Mitchell et al. 2014, Chiodi et al. 2018. Future climate scenarios for the Southeastern Coastal Plain include more frequent extreme weather events, such as hurricanes and severe droughts, as well as differences in the amount and timing of precipitation (IPCC 2014). Longleaf pine ecosystems have been hypothesized to exhibit higher resistance to these projected changes, leading to renewed efforts for their restoration. The success of these efforts will bring to the forefront the importance of accurate and unbiased estimates of their NEP.
Our results demonstrate large differences within savanna ecosystems in the contribution of the understory to ecosystem productivity and recovery, highlighting the critical need to evaluate how variation in savanna structure affects their contribution to global estimates of NPP. Furthermore, the inclusion of site-specific phenology would improve predictions regarding physiological patterns and ecosystem carbon dynamics, as represented in global dynamic vegetation models (Boke-Ol en et al. 2016). Further studies are needed to quantify these differences to fully understand how multiple disturbances interact over a range of structural variations in savanna ecosystems across the globe.

CONCLUSIONS
Understory phenology and fire are critical drivers of productivity in many ecosystems globally. We demonstrated that understory contribution to ecosystem fluxes is not uniform across environmental gradients within a savanna ecosystem, which affected ecosystem recovery times and overall productivity at the sites. As wildfire activities increase under changing climate (Westerling et al. 2006) and ecosystem management efforts increasingly focus on savanna ecosystems (Brudvig 2011, Noss et al. 2015, understanding climate and fire controls on understory carbon assimilation across edaphic variation will be critical for understanding ecosystem resilience.
For the longleaf pine savannas in this study, xeric sites characterized by lower overstory density were shown to have higher contributions of the understory to overall ecosystem fluxes, which increased their carbon sink capacity when understory NDVI was high. Native hardwood components to the understory of xeric sites exerted strong influence of phenological patterns observed in this study. More mesic and more productive sites experienced lower variability in ecosystem fluxes, as a result of the overstory dominating ecosystem fluxes, due to high basal and leaf area.
Site recovery from fire disturbance was context specific to interannual rainfall. For example, understories at more mesic sites recovered more rapidly from prescribed fires, which increased their carbon sink potential. In contrast, xeric sites showed little changes in NDVI and therefore NEP, GEP, and R eco following the fires, which resulted in little change in their sink capacity. In contrast, xeric site understories were less affected by drought, but did not improve overall ecosystem recovery from the disturbance.
Savanna ecosystems do not exhibit uniform ecosystem productivity and recovery from disturbance as a result of variations in structure from the underlying geology, which should be incorporated in global classifications of savanna ecosystems. We highlight the importance of incorporating site variations within savanna ecosystems to accurately predict ecosystem function and recovery from disturbances. Our results demonstrate the need for more fine-scale studies, especially in ecosystems with large structural variations, to accurately predict global carbon fluxes.