Carbon Dynamics of Pinus palustris Ecosystems Following Drought

Drought can affect forest structure and function at various spatial and temporal scales. Forest response and recovery from drought may be a result of position within landscape. Longleaf pine forests in the United States have been observed to reduce their carbon sequestration capacity during drought. We collected eddy covariance data at the ends of an edaphic longleaf pine gradient (xeric and mesic sites) over seven years; two years of normal rainfall were followed by 2.5 years of drought, then 2.5 years of normal or slightly above-average rainfall. Drought played a significant role in reducing the physiological capacity of the sites and was compounded when prescribed fire occurred during the same periods. The mesic site has a 40% greater basal area then the xeric site, which accounts for its larger sequestration capacity; however, both sites show the same range of variance in fluxes over the course of the study. Following drought, both sites became carbon sinks. However, the xeric site had a longer carry-over effect and never returned to pre-drought function. Although this study encompassed seven years, we argue that longer studies with greater spatial variance must be undertaken to develop a more comprehensive understanding of forest response to changing climate.


Introduction
Limitations in water availability in forested ecosystems can elicit varied responses over contrasting spatial and temporal scales.On short temporal and spatial scales, individual trees or stands can alter their leaf morphology to limit the loss of water from the system and try to reduce the stress that is being imposed on the plant and the ecosystem [1].Forests may also reduce their physiological activity, on the order of minutes to hours, to reduce water loss when demands are greatest during drought or high afternoon vapor pressure deficits [2].These reductions in physiological activity are typically associated with stomatal limitations, driven by lack of soil water availability.
Individual forests can respond to drought in a variety of ways.During seasonal droughts (several months in length), forests may facultatively reduce their leaf area with a premature dropping of foliage or a reduction in leaf development to reduce evaporative demands on the system [3][4][5].It has been hypothesized that the timing of leaf flush co-occurs with a plant's ability to maximize carbon Forests 2016, 7, 98 2 of 25 uptake prior to seasonal drought becoming so severe that respiration dominates the leaf-level carbon balance, which leads to an overall reduction in the plant's carbon stores [6].This strategy, ultimately, aids in the conservation of water and, at the same time, limits additional carbon losses; however, this also limits the plant's and forest's ability to sequester carbon immediately after the drought is alleviated [4,6].
Over long time periods and large spatial scales (100 to 1000 years and 100 to 1000 km 2 ), the landscape is a mosaic of smaller-scale ecosystems [7].These mosaics, in many cases, have developed in association with drought, which has also been shown to influence the structure and function of forested ecosystems based on their position within the landscape.The combination of drought, legacy effects, and underlying environmental conditions further drives gene expression, and leads to specific phylogenetic traits that create variations in forest types and their functions [7][8][9].This variation due to environmental pressures and landscape position has led to varying degrees of resilience in forest ecosystems along their distributional gradient [10,11].
With climate change, questions have arisen regarding the future distribution and survival of several key trees and forest types across the globe [12].Shifts in climate envelopes are expected to cause changes in tree species' distributional ranges [12][13][14], which will lead to greater uncertainty in estimates of future carbon sequestration of these forests.To better understand changing conditions and their roles on forest carbon capacity, the scientific community must undertake additional studies that utilize environmental gradients as "treatments" to facilitate the extrapolation of results to larger spatial scales and modeling efforts [15].Only then will we begin to develop the understanding needed to predict the future structure and function of forests [11,16,17].
We began a long-term study using an edaphic moisture gradient as a proxy for climate change in longleaf pine (Pinus palustris Mill.)-dominated forests within the Southeastern United States (SEUS), to better understand their carbon dynamics [18][19][20].The longleaf pine forest is an excellent example of how placement of a stand within the landscape has led to variation in structure and function, particularly when influenced by water limitations [21][22][23][24].Soil water holding capacity has been shown to be the primary variable that controls net primary production (NPP) within the longleaf ecosystem [21], and was linked to plant rooting zones and access to the water table [23,25].
Longleaf pine has historically been a significant part of the landscape in the SEUS, but intensive forest management has led to large-scale forest conversion to loblolly pine (Pinus taeda L.) plantation [26], resulting in a fragmented habitat.While the longleaf system will continue to be subjected to anthropogenic influence, this fire-maintained system is now dependent on humans to manage and maintain ecosystem integrity with low intensity, high frequency prescribed fire [27].Moreover, climate change is expected to alter temperatures and precipitation patterns in the SEUS region, which, in turn, may reduce managers' ability to maintain burn cycles, placing additional stress on the system and its ability to sequester carbon [14].
In this study, we follow up on the early findings of Whelan et al. [18,19], which examined the role that drought plays in the biogeochemical cycles of longleaf ecosystems along an edaphic gradient.These earlier studies showed that drought played a larger role in controlling the carbon dynamics of these forests compared to prescribed fire.Fire's effect on these systems was limited to ~30 days post-fire [20], after which, carbon uptake retuned to pre-fire rates.However, our earlier study also emphasized the need for longer-term studies that incorporate time series covering larger variations in climate.In our current study, we focus on the two sites at the ends of this edaphic moisture gradient, and follow their recovery from drought (as discussed in References [18,19]).This synthesis activity presents an opportunity to reduce uncertainties in our understanding of the pace, pattern, and process of carbon dynamics of differently adapted longleaf sites over a range of climates.Dependent on their position in the landscape, forest sites may respond differently to cumulative, chronic environmental stresses and disturbances [28].This study also sheds light on differences in adaptive strategies of forest ecosystems and their recovery capabilities as a result of their placement in the landscape.We hypothesize that (i) the magnitude of physiological change during recovery from drought would Forests 2016, 7, 98 3 of 25 be minimal in the xeric site because the plants' phenotypical responses have evolved with water as a limited resource, as compared to the mesic site, and (ii) following a drought, the mesic site would be a greater sink of carbon compared to the xeric site.The increase in the ecosystem's capacity to sequester more carbon would be associated with an increase in photosynthetic capacity following recovery from drought.
In this study, we tested the physiological variations of whole ecosystem processes (carbon capture) at each site under drought and non-drought conditions, with emphasis on drought recovery by examining and quantifying changes in light and temperature response on carbon processes over a seven-year period, which included ~4 years of 'normal' conditions and ~3 years of drought.We also utilized orthonormal wavelet transformation to analyze the frequency distribution of net ecosystem exchange of CO 2 (NEE), gross ecosystem CO 2 exchange (GEE) and ecosystem respiration (R eco ) rates during pre-drought, drought and recovery from drought.In addition, we examined diurnal activity of carbon dynamics to detect the physiological adaptive capacity of these two sites.

Study Sites
The study was conducted at the Joseph Jones Ecological Research Center (JJERC) in Southwestern Georgia, USA, from 22 October 2008 through 21 October 2015.The JJERC is an 11,000 ha reserve located within the Plains and Wiregrass Plains subsections of the Lower Coastal Plain and Flatwoods section.The center has extensive stands of second-growth longleaf pine, and has been managed with low intensity, dormant-season prescribed fires for the past ~75 years [21].The stands used for this study are approximately 90 years in age and were naturally regenerated.These stands have maintained their natural structure since they have been managed with single tree selection harvests using ecological forestry practices [29].
We selected study sites with soils in two drainage classes falling on the ends of an edaphic gradient for longleaf pine-wiregrass ecosystems at the center: (i) moderately poorly drained; and (ii) excessively well drained [30].The poorly drained sites occur on upland terraces with soils classified as Aquic Arenic Kandiudults.These soils are sandy loam over sandy clay loam or clay on nearly level slopes with a water holding capacity of 40 cm water m ´1 soil [21].An argillic horizon is present within 50 cm of the soil surface.The excessively well-drained sites occur on upland sand ridges of undulating slopes of 3%-4% and have deep, sandy soils, with no argillic horizon, i.e., no significant accumulation of clay is found within 300 cm of the surface.These soils are Typic Quartzipsamments with a water-holding capacity (in the upper 300 cm) of ~18 cm water m ´1 of soil [21].Hereafter, we refer to the moderately poorly drained site as mesic and excessively well-drained site as xeric.
The mesic ecosystem is characterized by a single overstory of dominant longleaf pine.Xeric ecosystems are dominated by open stands of longleaf pine and scrub oaks, where the co-dominant oak, turkey oak (Quercus laevis Walter), occurs in all the aboveground strata.Dense ground cover at all sites is dominated by the perennial grass, wiregrass (Aristida stricta Michx.), with numerous species of other perennial grasses and forbs also present [31].The sites used in this study have previously been described by the authors of [21,23,30,32].
Prior to the study, both sites were burned in winter 2007.During this study, both sites were burned in January 2009.The sites were on a two-year burn cycle and were both burned in late winter of 2011 and 2013.In 2015, the mesic site was burned in late winter, whereas the xeric site was burned in early spring.Each fire was set using backing fires on the downwind side of the management unit to move the fire away from the downwind fire-break.Strip head fires were then set perpendicular to the wind, starting on the downwind end of the unit and moving upwind.Fires were initiated ~30-50 m apart depending on fuel patterns, local weather conditions, and fire behavior [33].We limit our discussion of fire in this manuscript due to the short-duration, limited response these systems have shown to fire [20].[34,35].NEE is commonly estimated through simplification of the continuity equation by applying a control volume approach, where the integrals in Equation ( 1) are the vertical rate of change of mean molar CO 2 concentration (term I) and the vertical scalar flux divergence (term II) from ground level [36],

NEE was measured using open-path eddy covariance (EC) methods
where, the CO 2 concentration (C, µmol CO 2 m ´3) and the vertical velocities (w, m¨s ´1) were measured at a fixed plane above the plant canopy.In Equation ( 1), term II was estimated by the covariance of the turbulent fluctuations of C and w, where the turbulent fluctuations are the instantaneous deviation (at 10 Hz) from the mean block average (over 30 min).Primes denote the fluctuations; bars denote the mean from the 30-min averaging period.Micrometeorological convention is used here, where negative NEE values indicate ecosystem uptake of carbon.EC instruments were installed in early summer of 2008.
The EC measurement systems consisted of a three-dimensional sonic anemometer (CSAT-3, Campbell Scientific Inc., Logan, UT) and an open-path infrared gas analyzer (IRGA; LI-7500, LI-COR Biosciences, Lincoln, NE, USA) mounted at ~34.4 m (mesic) and ~34.9 m (xeric) at ground level (a.g.l.) on an antennae style tower.To minimize flow distortion between sensors, the IRGA and sonic anemometer were placed ~0.2 m apart, such that the open-optical path of the IRGA was vertically aligned to match the sonic volume height of the CSAT.Digital signals from the sonic anemometer and the gas analyzer (with factory 230 ms delay) were collected by a datalogger (CR3000, Campbell Scientific Inc., Logan, UT, USA) at 10 Hz.The IRGAs were calibrated every ~30 days using a traceable standard for CO 2 (600 ppm ˘1.0%), a dew point generator for H 2 O (LI-610, LI-COR Biosciences, Lincoln, NE, USA), and N 2 gas scrubbed with soda lime and Drierite, according to the protocols outlined by AmeriFlux [36].
Other soil meteorological measurements were collected on a second datalogger (model CR10X, Campbell Scientific Inc., Logan, UT, USA) and included soil heat flux plates buried at 8 cm (HFP01, Hukseflux Thermal Sensors, Delft, The Netherlands), volumetric water content (VWC) integrated over the top 30 cm of soil (CS616 water content reflectometer, Campbell Scientific, Logan, UT, USA), and soil temperature at 4 and 8 cm below the soil surface (Type-T copper/constantan thermocouples, Omega Engineering, Inc., Stamford, CT, USA).The aforementioned variables were measured every 10 s and averaged every 30 min.We also acquired Palmer Drought Severity Indices (PDSI) for the southwest Georgia region from the National Oceanic and Atmospheric Administration data archive [37].

Data Processing
Raw EC data were processed using EdiRe (version 1.4.3.1184;[38]), which carried out a 2-day coordinate rotation of the horizontal wind velocities to obtain turbulence statistics perpendicular to the local streamline.The covariance between turbulence and scalar concentrations was maximized through the examination of the time series at 0.1-s intervals on both sides of a fixed lag-time (in this case, ~0.3 s).Because of the relatively short roughness lengths and uniform canopy structure at these sites, we assumed that the influence of coherent structures and low frequency effects were captured by this approach.Fluxes were calculated for half-hour intervals and then corrected for the mass transfer resulting from changes in density not accounted for by the IRGA [39,40].Barometric pressure data were used to correct the fluxes to 1 standard atmosphere.
Flux data screening was applied to eliminate 30-min fluxes resulting from systematic errors, such as: (i) rain and condensation in the sampling path; (ii) incomplete 30-min datasets during system calibration or maintenance; (iii) poor coupling of the canopy with the external atmospheric conditions, as defined by the friction velocity, frictional velocity of wind (u*), using a threshold <0.20 m¨s ´1 [41,42]; and (iv) excessive variation from the half-hourly mean based on an analysis of standard deviations for u, v, and w wind and CO 2 statistics.Quality assurance of the flux data was also maintained by examining plausibility tests (e.g., ´30 < NEE < ´30 µmol m ´2¨s ´1), stationarity criteria, and integral turbulent statistics [43,44].
Eddy covariance measurements of NEE were estimated at a temporal resolution of one hour or less [45,46], such that: where, NEP is net ecosystem production and GPP is gross primary production.GPP cannot be measured directly, but rather is estimated from the right hand terms in Equation (2b).Half hourly fluxes of NEE (µmol m ´2¨s ´1) were used to calculate GEE and R eco in g C m ´2¨s ´1 from Equations (2a) and (2b) based on the logic provided in References [46][47][48].
To estimate seasonal and annual NEE, GEE, and R eco values, missing half hourly data were gap-filled using separate functions for day and night (NEE day , NEE night ).When photosynthetically active radiation (PAR) was ě10 µmol m ´2¨s ´1, daytime NEE data were gap-filled using a Michaelis-Menten approach [49]: and, when PAR was <10 µmol m ´2¨s ´1, nighttime NEE data were gap-filled using a modification of Lloyd and Taylor [50] approach: where, α is the apparent quantum efficiency (´µmol CO 2 µmol quanta ´1), ϕ is PAR (µmol quanta), R d is ecosystem respiration (µmol CO 2 m ´2¨s ´1), P max is the maximum ecosystem CO 2 uptake rate (µmol CO 2 m ´2¨s ´1), R 0 is the base respiration rate when air temperature is 0 ˝C, and b is an empirical coefficient (i.e., Van't Hoff's exponential function [51]).These functional relationships were estimated on a monthly basis to gap-fill where enough data were available.When too few observations were available to produce stable and biologically reasonable parameter estimates, equations estimated from annual data were used to gap-fill data by site.Gap-filled data accounted for 31% and 33% of daytime and 61% and 66% of nighttime values for mesic and xeric sites, respectively.Error estimation from gap-filled values of NEE was performed via bootstrap methods following Whelan et al. [18].Bootstrap procedures were performed monthly for gap-filling models to estimate missing NEE day data for all but two of the 168 site-months (the first two months at the xeric site), and for 135 of the 168 site-months to estimate missing NEE night data.All other months had too few usable Forests 2016, 7, 98 6 of 25 observations to estimate monthly equations.Therefore, we generated bootstrap samples for 15 annual equations (one daytime and 14 nighttime).In all cases, parameter estimates from the original data were within the 95% confidence region constructed from the bootstrap samples (Tables A1 and A2 Appendix A).

Data Analyses
Statistical models were formulated to test hypotheses about relationships between recovery from drought and CO 2 exchange, as well as the environmental drivers of those relationships.To evaluate the differences in physiological activity following before, during and after drought, we examined changes in the parameterization of NEE light response and temperature response curves over the time period.We also examined half-hourly physiological activity for each year to determine potential hysteresis.In order to quantify the underlying frequencies in NEE and R eco , we utilized orthonormal wavelet transformation (OWT), which is appropriate for time series such as these, since they do not require stationarity and allow for missing observations.We also employed OWT to analyze the frequency distributions of synchronously-measured NEE and R eco with micrometeorological drivers.
Potential hysteresis was graphically examined by averaging NEE at each half-hour at each site during each year of the study, from early morning to afternoon.Physiological activity was measured as the area within the polygon via the polygon function in ImageJ [52].Using a set value of pixels per cm and identical image sizes, the resulting measurements can be directly compared as relative values of light response hysteresis at the two sites for the years prior, during and post drought.
Response curves (Equations ( 3) and ( 4)) were fit over the entire time period of the study utilizing a moving window of 30 days to characterize changes in the response curves before, during, and after drought.The five parameters α, P max , R d , R 0 , and b were estimated via the SAS procedure PROC NLIN (SAS software, Version 9.2; SAS Institute Inc., Cary, NC, USA), exclusively with non-gap-filled observations.Estimated parameters were then used as response variables in five separate general linear models (GLMs) via the SAS procedure PROC MIXED.In order to isolate changes in the shape and asymptotic behavior of the curve due to drought response, model effects included PDSI for Southwest Georgia and the site.A covariate was also included to indicate if the 30-day moving window included prescribed fire within 30 days.This fire recovery period was employed, as previous results have shown the effects of fire are minimal on carbon dynamics after this point [20].In addition, a covariate for average temperature during the 30-day period was included as a proxy for seasonal effects.To meet the assumption of normality and homoscedasticity, each response variable was log-transformed prior to analysis.Due to the high temporal autocorrelation of the response variables, the assumption of independent observations was clearly violated, and therefore this analysis is presented in a descriptive context only.Marginal mean estimates of the curve parameters were used to graphically demonstrate the effect of site and drought on the NEE light and temperature response.
We employed OWT to NEE, GEE and R eco measurements following Reference [53].The authors of [53,54] present a conceptual description of the wavelet technique and application in flux research.(A detailed discussion of wavelet analysis and application in flux research can also be found in References [53,55,56]).OWT has been shown to be appropriate for spectral analyses of flux data, and is efficient in decomposing flux signals into both frequency and time-scale [53,56,57].Here, we use OWT to analyze the spectra of NEE, GEE and R eco during each of the 7 years of our study (2 pre-drought, 3 drought and 2 post-drought years).
OWT decomposes the given time series into two coefficient components (wavelet and scale) based on the finite basis function used [57].For our analysis, we computed orthonormal wavelet coefficients using the square-shaped Haar basis, due to its locality in the temporal domain and minimum effects of gaps or discontinuities [54].
To compute the OWT coefficients for NEE, we used strictly non-missing half-hourly observations.For GEE and R eco , all imputed values were used.Since OWT coefficients are computed at intervals of 2 n observations, and there are 17,520 = 2 14.097 half-hours in an annual dataset, annually-occurring cycles in 30-min flux data can be well-described with half-hourly data.However, diurnal cycles are not well-described since 48 = 2 5.585 .Therefore, we summarized half-hourly observations to 1.5 h, enabling computation of coefficients at time scales of 12 h, 24 h and ~1 month (32 days), and began each 24-h period for analyses at 6 a.m.Eastern Standard Time.The time series were then normalized to have zero mean and unit variance.Each year of corresponding GEE, NEE and R eco time series then underwent zero padding to fill in missing observations and to ensure that the number of data points was a power of 2, as the transformation decomposes the signal to half of the number of coefficients at each subsequent scale [55].The time series were then normalized again to have unit variance.The wavelet transformation was performed using 'wavelets' package in RStudio [58,59].OWT coefficients at each time-scale were summarized as the mean of squared coefficients to show the variance of the flux signal at each respective time-scale following Reference [55].

Micrometerological Conditions
Precipitation varied among sites and between years, with ~3 years having precipitation below the long term average ( [60], Table 1).During the study years of 2008-2009, 2012-2013 and 2013-2014, precipitation was near the long-term average for the area.The drought began in August of 2010 and continued until January of 2013.This period of drought led to PSDI values as low as ´5.2 (Figure 1).During this period, drought could be classified as moderate to extreme drought [61].Precipitation was distributed over the years evenly with a decrease in overall precipitation during the drought.There was one anomalous weather event in February 2013 with approximately ~20% of the annual rainfall occurring in one large rain event (Figure 2).Temperatures during the study were within the average range for the study area, with winter and summer average temperatures of 10.4 ˝C and 26.8 ˝C, respectively [60].

Annual Carbon
During years of severe drought, we observed declines in all three components of the carbon cycle (NEE, GEE and R eco ) with the largest declines in the 2010-2011 year, which corresponded to the lowest annual precipitation rates and PDSI values of <´5.0 (Figure 2).Interestingly, both sites seemed to have similar variance for their fluxes over the 7 years of the study; however, mean carbon uptake from the mesic site was large and the site remained a sink over the entire study period (Table 1).During periods of drought the xeric site carbon balance was driven by both a reduction in gross photosynthetic rates and by a dominance of respiration rates over NEE day (Table 1).These results are further validated when examining the changes in averaged diurnal patterns of NEE, GEE and R eco during recovery from drought when physiological activity increased at both sites (Figure 1).

Annual Carbon
During years of severe drought, we observed declines in all three components of the carbon cycle (NEE, GEE and Reco) with the largest declines in the 2010-2011 year, which corresponded to the lowest annual precipitation rates and PDSI values of <−5.0 (Figure 2).Interestingly, both sites seemed to have similar variance for their fluxes over the 7 years of the study; however, mean carbon uptake from the mesic site was large and the site remained a sink over the entire study period (Table 1).During periods of drought the xeric site carbon balance was driven by both a reduction in gross photosynthetic rates and by a dominance of respiration rates over NEEday (Table 1).These results are further validated when examining the changes in averaged diurnal patterns of NEE, GEE and Reco during recovery from drought when physiological activity increased at both sites (Figure 1).

Light and Temperature Responses Following Drought
When adjusted for the seasonality of temperature, the apparent quantum efficiency, α, responded in opposite ways to drought (as measured by PDSI) when comparing both the xeric versus mesic site, and the fire recovery versus post-recovery period (Figure 3A).During the fire

Light and Temperature Responses Following Drought
When adjusted for the seasonality of temperature, the apparent quantum efficiency, α, responded in opposite ways to drought (as measured by PDSI) when comparing both the xeric versus mesic site, and the fire recovery versus post-recovery period (Figure 3A).During the fire recovery period (i.e., observations within 30-day post-fire), drier conditions prevailed (and decreasing PDSI, moving right to left in Figure 3A), and were associated with less negative values of α at the xeric site but more negative α at the mesic site, resulting in higher capacity for C uptake (reference Equation (3) and micrometeorological convention).Conversely, during the post-recovery period (i.e., observations outside of the 30-days post-fire), drier conditions were associated with moderately lower α in the xeric site, and the relationship between PDSI and α closely resembled that of the mesic site during the fire recovery period, while at the mesic site, drier conditions were associated with slightly higher α.
The maximum ecosystem CO 2 uptake rate, P max , responded similarly to changes in PDSI when adjusted for the seasonality of temperature (Figure 3B).P max was higher under drier conditions (low values of PDSI) at both sites and with all fire conditions.However, the response of P max to drier conditions was markedly stronger in the fire recovery period at the xeric site.
There was very little change in the ecosystem respiration parameter, R d , associated with changes in PDSI during the fire recovery period (Figure 3C).During the fire recovery time, R d responded to PDSI in opposite ways depending on site.R d was lower during drier conditions at the xeric site, but higher during drier conditions at the mesic site.
Forests 2016, 7, 98 9 of 27 (reference Equation ( 3) and micrometeorological convention).Conversely, during the post-recovery period (i.e., observations outside of the 30-days post-fire), drier conditions were associated with moderately lower α in the xeric site, and the relationship between PDSI and α closely resembled that of the mesic site during the fire recovery period, while at the mesic site, drier conditions were associated with slightly higher α.
The maximum ecosystem CO2 uptake rate, Pmax, responded similarly to changes in PDSI when adjusted for the seasonality of temperature (Figure 3B).Pmax was higher under drier conditions (low values of PDSI) at both sites and with all fire conditions.However, the response of Pmax to drier conditions was markedly stronger in the fire recovery period at the xeric site.
There was very little change in the ecosystem respiration parameter, Rd, associated with changes in PDSI during the fire recovery period (Figure 3C).During the fire recovery time, Rd responded to PDSI in opposite ways depending on site.Rd was lower during drier conditions at the xeric site, but higher during drier conditions at the mesic site.When adjusted for the seasonality of temperature, the response of R0, the base respiration rate when air temperature is 0 °C, to changes in PDSI was similar to that of Rd during the fire recovery period (Figure 4A), with slightly higher values of R0 during drier conditions.During the fire recovery period, R0 responded similarly to drier conditions at the mesic site, but was slightly lower When adjusted for the seasonality of temperature, the response of R 0 , the base respiration rate when air temperature is 0 ˝C, to changes in PDSI was similar to that of R d during the fire recovery period (Figure 4A), with slightly higher values of R 0 during drier conditions.During the fire recovery period, R 0 responded similarly to drier conditions at the mesic site, but was slightly lower in magnitude when compared to the post-fire recovery period.At the xeric site, R 0 had a much stronger response to PDSI during the fire recovery period; R 0 was much lower during drier conditions.
At the xeric site, the empirical coefficient b also had a strong response to PDSI during the fire recovery period; b was much higher during drier conditions (Figure 4B).In contrast, b was slightly lower during drier conditions at the mesic site.During the post-recovery period, b was also lower with drier conditions. in magnitude when compared to the post-fire recovery period.At the xeric site, R0 had a much stronger response to PDSI during the fire recovery period; R0 was much lower during drier conditions.
At the xeric site, the empirical coefficient b also had a strong response to PDSI during the fire recovery period; b was much higher during drier conditions (Figure 4B).In contrast, b was slightly lower during drier conditions at the mesic site.During the post-recovery period, b was also lower with drier conditions.

Orthonormal Wavelet Analyses
At the mesic site, NEE spectra for each year showed that variability at the 12-h time scale was higher during pre-drought years versus post-drought years, while higher variability was observed during all three drought years and the second post-drought year at 24-h time scale (Figure 5A).The change in spectral peaks from the pre-drought years to the drought years also occurred progressively, as evident by a transition of wider spectral peak during the second pre-drought year (at both 12-h and 24-h time-scales).This mesic pattern was mostly matched at the xeric site (Figure 5B).However, at the bi-weekly to monthly time-scales, as well as subsequent time-scales during the second drought year, the xeric site had higher variability compared to the mesic site.In addition, during the first post-drought year, variability at 12-h-24-h was lower at the xeric site than at the mesic site.At both sites, the spectra of the second post-drought years follow the same pattern as those of the other drought year spectra at higher frequencies at most time-scales (Figure 5A,B).

Orthonormal Wavelet Analyses
At the mesic site, NEE spectra for each year showed that variability at the 12-h time scale was higher during pre-drought years versus post-drought years, while higher variability was observed during all three drought years and the second post-drought year at 24-h time scale (Figure 5A).The change in spectral peaks from the pre-drought years to the drought years also occurred progressively, as evident by a transition of wider spectral peak during the second pre-drought year (at both 12-h and 24-h time-scales).This mesic pattern was mostly matched at the xeric site (Figure 5B).However, at the bi-weekly to monthly time-scales, as well as subsequent time-scales during the second drought year, the xeric site had higher variability compared to the mesic site.In addition, during the first post-drought year, variability at 12-h-24-h was lower at the xeric site than at the mesic site.At both sites, the spectra of the second post-drought years follow the same pattern as those of the other drought year spectra at higher frequencies at most time-scales (Figure 5A,B).
The most striking observation for GEE was its lower magnitude spectrum for the first post-drought year at both sites.The spectral gap (weekly to monthly time-scales; as described by Stoy et al. [53]) was comparatively lower than other years (Figure 5C,D).In the case of R eco , the spectra were almost identical among sites except during the third drought year, where the mesic site had a higher variance (Figure 5E,F).
The most striking observation for GEE was its lower magnitude spectrum for the first post-drought year at both sites.The spectral gap (weekly to monthly time-scales; as described by Stoy et al. [53]) was comparatively lower than other years (Figure 5C,D).In the case of Reco, the spectra were almost identical among sites except during the third drought year, where the mesic site had a higher variance (Figure 5E,F).

NEE Activity
The NEE activity period shifted progressively from ~7:30-18:30 to ~6:30-17:30 during drought years at both sites.However, the recovery from drought to post-drought was stronger (with a higher magnitude of maximum photosynthetic assimilation) at the xeric site than at the mesic site (Figure 6).The NEE in the first post-drought year at the xeric site showed proportionally (to pre-and post-drought years) higher magnitude of maximum photosynthetic assimilation, while that of the mesic site showed proportionally (to pre-drought years) lower magnitude during the first drought year (Figure 6).The relative integrated area increased at both sites until the second drought year, indicating reductions in daytime physiological activity.During the third drought year and post-drought period, the area at the mesic site decreased, indicating an increase in the length of daytime activity, while at the xeric site, activity remained at the reduced level until the second post-drought year.Reco showed higher magnitude during the second post-drought year than during other years at the xeric site (Figure 6E,F).The mesic site post-drought years had higher magnitudes

NEE Activity
The NEE activity period shifted progressively from ~7:30-18:30 to ~6:30-17:30 during drought years at both sites.However, the recovery from drought to post-drought was stronger (with a higher magnitude of maximum photosynthetic assimilation) at the xeric site than at the mesic site (Figure 6).The NEE in the first post-drought year at the xeric site showed proportionally (to pre-and post-drought years) higher magnitude of maximum photosynthetic assimilation, while that of the mesic site showed proportionally (to pre-drought years) lower magnitude during the first drought year (Figure 6).The relative integrated area increased at both sites until the second drought year, indicating reductions in daytime physiological activity.During the third drought year and post-drought period, the area at the mesic site decreased, indicating an increase in the length of daytime activity, while at the xeric site, activity remained at the reduced level until the second post-drought year.R eco showed higher magnitude during the second post-drought year than during other years at the xeric site (Figure 6E,F).The mesic site post-drought years had higher magnitudes of maximum R eco than during the second pre-drought year, but less than that of the first pre-drought year. of maximum Reco than during the second pre-drought year, but less than that of the first pre-drought year.

Discussion
P. Palustris has a large geographic range in the Southeastern US that also represents a large environmental gradient.In this study, two sources of disturbance, drought and fire, and subsequent recovery were superimposed on this environmental gradient.Large ranges in ecosystem functions were observed among our sites, and through the disturbances and recovery; yet we did not observe any increased mortality rates [18].We attribute this large range of plasticity in ecosystem functions and ability to recover to phenotypic expression of the large gene flow [62].Woody plant gymnosperms, wind-pollinated woody species and species with larger geographical ranges have much higher gene flow than their angiosperm counterparts [63].We do not attribute this functional plasticity to natural selection, but cannot rule it out, particularly in light of future and chronic disturbance of climate change [28].

Disturbance and Recovery
P. Palustris ecosystems have co-evolved with fire and drought.Natural-and now managed-high frequency and low intensity burns have been shown to support the high degree of biodiversity and endemism, and support the intermediate disturbance hypothesis [32].Hence, one would expect high degrees of biodiversity and associated diversity of plant functions would scale to the ecosystem level and contribute towards resiliency and the ability to recover after disturbance [64,65].Indeed, Whelan et al. [19] show that these ecosystems return post-fire to pre-fire rates of net C uptake at ~30 days.Even still, these ecosystems require the ability to function and thrive with both fire and drought, even when the discrete timescales of these disturbances differ; with fire and drought being relatively shorter and longer, respectively.Moreover, fire and drought may occur consecutively in time, but more often than not, they occur concurrently.We recognize that the responses to fire are superimposed on that of drought, but because of the shorter recovery time to fire (i.e., 30 days), recovery of fire appear to be discounted from the longer term responses to drought.

Fire
As previously discussed in Whelan et al. [18,19] longleaf ecosystems recover from fire within a short time period (Figure 3A).This recovery has been thought to be associated with these ecosystems' ability to allocate a larger portion of their carbon reserves below ground [66].In addition, the low intensity fire limits crown damage lending to the chief component of NEE to remain intact, and limiting the declines in carbon uptake for the ecosystem [27].Fire is also introduced into these systems during the known growing season when understory grasses are dormant; the regrowth of grasses begins following fire and contributes to the systems rapid recover.The combination of these factors aids in the resiliency of the ecosystem to short-term perturbations such as fire, and buffers the ecosystem capacity for sequestering carbon.In turn, this has led to the discussion of the role that longer term disturbances, such as drought, play in the carbon sequestration capacity of these systems.

Drought
Using measures of productivity as an indication of resiliency [67], recovery from drought was dependent on the temporal scale in which we assessed our data.On the 24-h period, the maximum diurnal uptake rates were slightly less from the xeric site as compared to the mesic site for all drought and non-drought periods (Figures 1 and 6).We primarily attribute this result to the fact that the xeric site had ~40% less basal area.However, the diurnal patterns also differed, with a sharper decline in NEE rates during the afternoon at the xeric sites, i.e., hysteresis (Figures 1 and 6), resulting in less net C uptake in the xeric site when integrating under the daytime hours (Figure 6).We attribute this hysteresis in both sites to stomata closure taking place in response to increasing mid-day VPD (VPD threshold of 2.5 KPa detected by Whelan et al. [19]) and differences in soil water availability, together suggesting these trees still have access to deep water.P. Palustris's ability to change its C uptake in response to drought (and its recovery) is through changes in needle photosynthetic capacity (2-3 years needle cohort) and whole tree conductance [68].The reduction in NEE rates was more pronounced in the first year of the drought compared to the second year, even though the PDSI was the same or greater in the second year of drought (Table 1 and Figure 2).The similar pattern in NEE rates was observed in the recovery from drought, where the rates from the first year after drought were not as large as those from the second year.
Interestingly, the spectral patterns of NEE differed from pre-drought to drought conditions, where the amount of variance changed temporally (i.e., spectral density) shifting from a 12-h to a 24-h spectral peak (Figure 5).The diurnal R eco spectral peak was not altered pre-, post-or during the drought, and was ~12 h (Figure 5).This suggests that greatest amount of (diurnal) variance in R eco occurred every 12-h, as would be expected for nighttime respiration.However, the shift from 12-h peak spectral density in NEE suggests that the maximum (diurnal) variance changed from controls on nighttime and daytime processes, to the importance of day-to-day variability.Put another way, the day-to-day variance in NEE during drought periods-and the first year of recovery from drought-exceeded that of the variability observed in nighttime R eco and that observed in NEE day .However, because the R eco spectra did not change we attribute the shift in NEE spectra to that of NEE day alone.It also appears (Figure 5) that during year 2 post-drought, the NEE spectra are returning towards the 12-h peak.Taken in concert, these changes in patterns also speak to the annual phenotypic plasticity of this ecosystem.
There was a large range in annual productivity among the pre-, post and drought years.Even with prescribed burns and a two-year severe drought, these ecosystems were still a sink of C with mean annual NEE rates of 2.1 ˘0.18 and 0.7 ˘0.18 Mt C ha ´1¨year ´1 (µ ˘1 SE) for mesic and xeric sites, respectively, and a total C accumulation of 14.6 and 5.2 Mt C ha ´1 over the seven years of this study for the mesic and xeric sites, respectively.While the magnitude of GEE, NEE, R eco was larger at the mesic site as compared to the xeric, there was no significant difference when taking into account their variance.Hence, we could not reject Hypothesis 1: The magnitude of physiological change during recovery from drought would be minimal at the xeric site because the plants' phenotypical response has evolved with water as a limited resource, as compared to the mesic site.However, we could accept Hypothesis 2: Following a drought, the mesic site would be a greater sink of carbon compared to the xeric site.However, NEE is a very small difference between two large processes-GEE and R eco , and GEE cannot be measured directly.Hence, any discussion on the controls of GEE is then circular.
In 2009-2011, both sites became near C neutral, and our data actually suggests that the xeric site became a source of carbon.Annual R eco rates were of different magnitude among the mesic and xeric sites, but both exhibited similar variance (Table 1, Figures 1 and 2).This suggests that the annual controls on R eco are similar among sites.But because R eco is made up of both auto-and heterotrophic processes and the available C pools by which they rely, what could be the possible mechanism for R eco to exceed GEE?
Temporal scales aside, ecosystem respiration has been shown to be a function of temperature, and also available water in the case of below-ground respiration.Because there were small increases in aboveground biomass and seasonal leaf area dynamics did not change among years, we discount any new and large contributions to aboveground respiration pools towards R eco .One can imagine that there can be a delta (change) in aboveground maintenance and growth respiration to address additional 'cost' in response to drought, e.g., a delta to repackage the photosynthetic machinery in the needle.However, we assume these delta costs would be small, and we do not know of any study that has addressed this issue.Hence, from year-to-year, we assume the aboveground respiration is a function of available substrate, i.e., flat tax argument and the largest aboveground year-to-year change affecting productivity was associated with changing the water use efficiency (gram water used per gram C captured) of the new cohort of needles.
The ability for this ecosystem to enhance the capture of C primarily lies (first order) with annual replacement of old cohort of needles with new ones that have a drought-adapted capacity.This is an annual response, and with a full turnover of capability occurring with the full replacement of 2-3 cohorts of needles.We observe this annual lag in NEE time series (Figure 2) and the reduction of 'carry over' effects in year 2 of the drought, and year 2 of the post drought recovery.
Patterns of productivity were similar in the mesic and xeric sites, and exhibited similar temporal variances, but with generally larger mean quantities found in the mesic site.In addition, while respiration persisted pre-, post-and during the drought, both ecosystems were a net sink over seven-year and we did not detect any changes to mortality rates [19].We show a large range in phenotypic plasticity in ecosystem functions which constrains a set of responses to disturbance, in this case drought.This large range of plasticity enables this ecosystem to continue to function, while retaining essentially the same function, structure, feedbacks, and identity, i.e., resiliency [69] and also has the capacity for the component ecosystem processes to manage resilience, i.e., adaptability [70].

Conclusions
IPCC [13] predicts that temperatures may increase up to 1.5 ˝C for the SE-US in combination with a 10% increase in precipitation and a decline in frost events for the region.However, the timing of rain events is also expect to be altered, with increases in heavy precipitation events and great periods of drought between these events.The combination of these changes leads to additional scientific questions regarding the future carbon cycles of longleaf pine forest and their ability to maintain structure and function [14].Our findings show that drought leads to stomatal limitations and legacy effects that alter carbon dynamics for at least one year following drought.Increases in drought condition in the future may further reduce the ecosystem's ability to sequester carbon.The increase in temperature may also increase the metabolic activity of the system leading to great water use [18,19].Changes in drought patterns may also alter management practices by reducing our ability to use fire as a prescription to maintain these ecosystems [14].The other consideration scientist needs to make is the possibility of carbon dioxide (CO 2 ) fertilization and how the combination of temperature and disturbance will alter a longleaf pines responses in the future [71].These alterations in environment, atmospheric CO 2 concentration.and fire regimes are expected to increase variability in a system's carbon sequestration capacity on an annual time scale [71].Our study, however, provides foundational evidence that longleaf forests, with their broad phenological plasticity have the ability to survive potentials changes in climate and management.On the other hand, these finding have been developed from one region of the distributional ranges for longleaf ecosystems.To reduce uncertainties within our understanding of the pace, patterns and process of the system, new studies must be undertaken.These studies must encompass a greater range of the distributional and environmental conditions in which longleaf pine forests exist if the scientific community is to gain a greater understanding of the ecosystem's survival with change.1.100 0.968 1.237 1.100 0.068 0.068 LCL = lower limit of 95% confidence region, UCL = upper limit of 95% confidence region.R 0 is the base respiration rate when air temperature is 0 ˝C and b is an empirical coefficient.

Figure 2 .
Figure 2. Monthly sums of Gross Ecosystem Exchange, Net Ecosystem Exchange and Ecosystem respiration for the mesic (A) and xeric sites (B), respectively; Palmer Drought Severity Index and summed monthly precipitation for the two stands (C).Shaded area indicates period of drought (PDSI <´3.0), and arrows indicate dates of prescribed fire occurrence.

Figure 2 .
Figure 2. Monthly sums of Gross Ecosystem Exchange, Net Ecosystem Exchange and Ecosystem respiration for the mesic (A) and xeric sites (B), respectively.Palmer Drought Severity Index and summed monthly precipitation for the two stands (C).Shaded area indicates period of drought (PDSI <−3.0), and arrows indicate dates of prescribed fire occurrence.

Figure 3 .
Figure 3. Least square mean estimated light response parameter values for fire recovery period and post-recovery by site: (A) Apparent quantum efficiency (α), (B) maximum ecosystem CO2 uptake rate (Pmax), and (C) ecosystem respiration (Rd).

Figure 3 .
Figure 3. Least square mean estimated light response parameter values for fire recovery period and post-recovery by site: (A) Apparent quantum efficiency (α), (B) maximum ecosystem CO 2 uptake rate (P max ), and (C) ecosystem respiration (R d ).

Figure 4 .
Figure 4. Least square mean estimated temperature response parameter values for fire recovery period and post-recovery by site (A) the base respiration rate when air temperature is 0 °C (R0), and (B) the empirical coefficient, b.

Figure 4 .
Figure 4. Least square mean estimated temperature response parameter values for fire recovery period and post-recovery by site (A) the base respiration rate when air temperature is 0 ˝C (R 0 ), and (B) the empirical coefficient, b.

Figure 5 .
Figure 5. Mesic versus xeric normalized wavelet spectra of net ecosystem exchange of CO 2 (NEE; (A,B)), gross ecosystem exchange (GEE; (C,D)), and ecosystem respiration (R eco ; (E,F)).The graphs show variance of the time-series signal at each respective time-scale.

Figure 6 .
Figure 6.Annual diurnal patterns of net ecosystem exchange of CO 2 (NEE) averaged at each half-hour in each site during each year of the study, from early morning to afternoon.Integrated area represents daytime physiological activity, relative to mesic site pre-drought year 1.

Table 1 .
Annual precipitation (mm), Net Ecosystem Exchange of Carbon (NEE, g C m 2 ¨year ´1), Gross Ecosystem Exchange of Carbon (GEE), Ecosystem Respiration (R eco ), Ratio of GEE to R eco .

Table A1 .
Cont. limit of 95% confidence region, UCL = upper limit of 95% confidence region.α is the apparent quantum efficiency (-µmol CO 2 µmol quanta ´1), R d is ecosystem respiration (µmol CO 2 m ´2¨s ´1), and P max is the maximum ecosystem CO 2 uptake rate.

Table A2 .
Distribution of parameters from nighttime net ecosystem exchange (NEE night ) bootstrap simulations.