El Niño Southern Oscillation (ENSO) Enhances CO2 Exchange Rates in Freshwater Marsh Ecosystems in the Florida Everglades

This research examines the relationships between El Niño Southern Oscillation (ENSO), water level, precipitation patterns and carbon dioxide (CO2) exchange rates in the freshwater wetland ecosystems of the Florida Everglades. Data was obtained over a 5-year study period (2009–2013) from two freshwater marsh sites located in Everglades National Park that differ in hydrology. At the short-hydroperiod site (Taylor Slough; TS) and the long-hydroperiod site (Shark River Slough; SRS) fluctuations in precipitation patterns occurred with changes in ENSO phase, suggesting that extreme ENSO phases alter Everglades hydrology which is known to have a substantial influence on ecosystem carbon dynamics. Variations in both ENSO phase and annual net CO2 exchange rates co-occurred with changes in wet and dry season length and intensity. Combined with site-specific seasonality in CO2 exchanges rates, El Niño and La Niña phases magnified season intensity and CO2 exchange rates at both sites. At TS, net CO2 uptake rates were higher in the dry season, whereas SRS had greater rates of carbon sequestration during the wet season. As La Niña phases were concurrent with drought years and extended dry seasons, TS became a greater sink for CO2 on an annual basis (−11 to −110 g CO2 m−2 yr−1) compared to El Niño and neutral years (−5 to −43.5 g CO2 m−2 yr−1). SRS was a small source for CO2 annually (1.81 to 80 g CO2 m−2 yr−1) except in one exceptionally wet year that was associated with an El Niño phase (−16 g CO2 m−2 yr−1). Considering that future climate predictions suggest a higher frequency and intensity in El Niño and La Niña phases, these results indicate that changes in extreme ENSO phases will significantly alter CO2 dynamics in the Florida Everglades.


Introduction
Teleconnections from the El Niño Southern Oscillation (ENSO) are known to strongly affect climate patterns across North America [1,2,3,4]. ENSO cycles are alternating periods of warm (El Niño phase) and cold (La Niña phase) Pacific Ocean surface temperatures [5,6], and have occurred with regular periodicity (3 to 7 years) over the last 130,000 years [3]. Shown to influence worldwide precipitation patterns [3], ENSO phases are also correlated with global terrestrial productivity [7] and climate anomalies [3,8].
In the Florida Everglades, changes in the long-term hydrologic cycle have been linked to extreme ENSO phases (El Niño and La Niña phases) [2,9]. Precipitation patterns in this region form wet and dry seasons, the frequency and magnitude of which fluctuate with changing climate patterns [8]. Here, El Niño phases increase dry season rainfall causing higher seasonal and annual water levels [2,9]. In contrast, La Niña phases reduce dry season rainfall, leading to extreme drought and the water table dropping below the soils surface [2,3,4,9]. Because annual shifts in carbon dioxide (CO 2 ) exchange rates have been linked to changes in surface hydrology in short-term studies [10,11,12,13], El Niño and La Niña phases may be an important driver of seasonal-to-interannual variations in hydrology and ultimately the productivity of Everglades freshwater marsh ecosystems. It is well known that wetland ecosystem structure and function is tightly coupled to hydrology, and as such it controls wetland carbon (C) sequestration [8,13,14,15]. Wetland CO 2 exchange rates respond to changes in surface hydrology [11,13,16]. The magnitudes of intra-and inter-annual fluctuations in surface hydrology are sensitive to global climate cycles [2], and directly affect CO 2 exchange. As a result, inter-and intra-annual fluctuations in CO 2 exchange rates in the Everglades region may be significantly influenced by El Niño and La Niña phases.
Increased atmospheric concentrations of CO 2 and other greenhouse gases are expected to alter the frequency of El Niño and La Niña phases [17]. In addition to the El Niño and La Niña-driven effects, climate change projections also suggest changes in the magnitude and frequency of seasonal precipitation patterns, as well as higher dry season temperatures [9,18]. Precipitation projections suggest wetter summers (wet season) and more severe drought (dry season) over the southeastern U.S. [19]. Fluctuations in water availability as a result of these changes may alter ecosystem structure and function.
Surface hydrology is managed differently among watersheds within Everglades National Park, which provides a unique opportunity to examine ecosystem function with differing hydroperiods, while still experiencing similar climate. Schedlbauer et al. [12] and Jimenez et al. [11] have assessed the effects of managed hydroperiods on seasonal and annual carbon (C) dynamics for short periods (1 to 2 years) in Everglades freshwater marsh ecosystems. However, there has been no research to date that has assessed the effects of ENSO teleconnections on seasonal and annual CO 2 dynamics in the Everglades. Because El Niño and La Niña phases are expected to alter the frequency and intensity of precipitation and temperature regimes, it is unknown how, when, and with what magnitude ecosystem CO 2 exchange rates will respond to these fluctuations. However, this information is key to develop a prognostic understanding of how these ecosystems will behave in the future.
The goal of this research is to understand the relationship between extreme ENSO phases and intra-and inter-annual fluctuations in CO 2 exchange rates (NEE, R eco , and GEE). We hypothesize that El Niño and La Niña will amplify the site-specific seasonal responses in CO 2 fluxes. At the short-hydroperiod site (Taylor Slough; TS) it has been shown that enhanced net carbon uptake rates are associated with the dry season, while at the long-hydroperiod site (Shark River Slough; SRS) greater net carbon uptake rates are associated with wet season conditions [11,16]. As El Niño and La Niña phases increase wet and dry season intensity in the Everglades region [2,3,4,9], we expect the site-specific seasonal response to change correspondingly with changes in season intensity. Season intensity here refers to deviations from the mean water availability, so that larger absolute numbers indicate the season intensity and the sign indicates wetter (+) or dryer (-) conditions. We also hypothesize that the variation in season length will explain differences in the interannual CO 2 dynamics in the respective processes of uptake and efflux. For example, longer and wetter wet seasons will increase the capacity to uptake carbon at long hydroperiod sites, whereas longer and hotter dry seasons will increase the capacity for carbon uptake at short hydroperiod sites. In this study we used the eddy covariance method to estimate whole ecosystem exchanges in CO 2 , and a combination of linear, non-linear and time series modeling techniques to statistically address these hypotheses.

Study Site
The Florida Everglades are classified as subtropical wetlands with a year-long growing season and distinct wet and dry seasons that define annual variation [20,21]. Water enters the Everglades through local precipitation events, which average 1380 mm annually [8], and through regional runoff. Presently, water dynamics are controlled by the South Florida Water Management District, which uses a complex system of canals, levees, and pumping stations [8,22]. The majority of rainfall (,70%) occurs during the wet season, which begins in May or early June with convective events and tropical depressions, e.g., thunderstorms and hurricanes [8]. Surface water levels generally increase throughout the wet season, are highest at wet season end in October, and decline to their lowest levels by dry season end in May [3]. During the dry season, the Bermuda High-pressure cell prevents convective clouds from forming thunderstorms, making continental fronts the main source of precipitation [20]. This switch from wet season tropical climate to dry season temperate climate causes distinct changes in the amount of precipitation in the region [20]. Dry season precipitation accounts for ,30% of annual precipitation [8].
The study sites are two oligotrophic freshwater marsh ecosystems that are within the Florida Coastal Everglades (FCE) long-term ecological research (LTER) program in Everglades National Park (FCE-LTER, http://fcelter.fiu.edu/research/ sites/; Fig. 1). Taylor Slough (25˚26916.50 N, 80˚35940.680 W) is a shorthydroperiod marsh that is inundated for 4 to 6 months each year (,June to November) and is characterized by shallow (,0.14 m) marl soils overlying limestone bedrock. Mean canopy height (Z) and surface roughness (d) for this site are 0.73 and ,0.3 m, respectively. Shark River Slough (25˚3396.720N, 80˚46957.360W) is a long-hydroperiod marsh that is inundated ,12 months each year and is characterized by peat soils (,1 m thick) overlying limestone bedrock with ridge and slough microtopography [23]. For this site, Z and d are 1.02 and ,0.4 m, respectively. Differences in hydroperiod result from spatial variability in elevation [3] and exposure to surface runoff.

Eddy Covariance and Meteorological Data
At each site, open-path infrared gas analyzers (IRGA, LI-7500, Li-COR Inc., Lincoln, NE) were used to measure CO 2 (c; mg mol -1 ) and water vapor molar density (r v ; mg mol -1 ), and a paired sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, UT) was employed to measure sonic temperature (T s ; K) and 3-dimensional wind speed (u, v and w, respectively; m s -1 ). These paired sensors were 0.09 m apart and installed at 3.30 and 3.24 m above ground level (a.g.l.) at TS and SRS, respectively. Data were logged at 10 Hz on a datalogger (CR1000, Campbell Scientific Inc.) and stored on 2 GB CompactFlash cards. Both IRGAs were calibrated monthly using a trace gas standard for CO 2 in air (+1.0%), dry N 2 gas and a portable dewpoint generator (LI-610, LI-COR Inc.). Footprint analyses [26,27] indicated that 80% of measured fluxes were within 100 m of the tower during convective conditions at both sites. Other meteorological variables were measured at 1-sec and collected as half-hourly averages, acquired by the same datalogger, and included: air temperature, (T air ;˚C) and relative humidity (Rh; %) (HMP45C, Vaisala, Helsinki, Finland) mounted within an aspirated shield (43502, R.M. Young Co., Traverse City, MI), and barometric pressure (P; atm) (PTB110, Vaisala). The T air /Rh sensors were installed at the same height a.g.l. as the IRGA and CSAT.
At each site, additional meteorological data was measured at 15-sec, and collected as 30-min averages through a multiplexer (AM16/32A Campbell Scientific Inc.) with another datalogger (CR10X Campbell Scientific Inc. Zonen). Precipitation measurements were made with tipping bucket rain gages (mm) (TE525, Texas Electronics Inc., Dallas, TX). Soil volumetric water content (VWC; %) was calculated from equations developed for peat and marl soils using the methodology of Veldkamp & O'Brien [28], from the dielectric constant using two soil moisture sensors (CS616, Campbell Scientific Inc.) installed at a 45˚angle at the soil surface, at each site. Soil temperature (T s ;˚C) was measured at 5 cm, 10 cm, and 20 cm depths at two locations within each site using insulated thermocouples (Type-T, Omega Engineering Inc., Stamford, CT). When inundated at SRS, water temperature, (T w ;˚C) was measured using two pairs of insulated thermocouples (Type-T, Omega Engineering Inc.), each pair located at a fixed height 5 cm above the soil surface and another attached to shielded floats that held the thermocouples 5 cm below the water surface. At TS, T w was measured using insulated thermocouples (Type-T, Omega Engineering Inc.) located at a fixed height 2 cm below the water surface. Water level (m) at both sites was recorded every half-hour with a water level logger (HOBO U20-001-01, Onset, Bourne, MA).

Data processing
Net ecosystem exchange (NEE) of CO 2 was estimated through simplification of the continuity equation by applying a control volume approach from the ground level to the top measurement height (z; m) [29]. Vertical windspeed (w) was first estimated mean to streamline using a 2-d rotation in a Cartesian coordinate framework [29]. NEE (mmol m 22 s 21 ) was then estimated using the covariance of the turbulent fluctuations of the vertical rate of change of mean molar density of CO 2 (c9), and the vertical scalar flux divergence (w9), where the turbulent fluctuations are the instantaneous deviation (at 10 Hz) from the mean block average (term I) over 30 min, and the storage flux (term II): where,w 0 c 0 is the measured covariance (m s 21 mmol C mol 21 ) of the molar density of CO 2 measured at a fixed plane above the plant canopy (Z), z is the vertical dimension, and the overbar is the averaging period, in this case is 30-min. NEE was then divided by the molar volume of air, V, (m 3 mol 21 ) to convert the units from density to molar fraction, i.e., mmol CO 2 m 22 s 21 , such that: where, R is the ideal gas constant (0.082 L atm K 21 mol 21 ), P is atmospheric pressure (1.10325 atm), and T k is the actual air temperature, estimated by: where, q is the molar fraction of water vapor calculated by unit conversion of r v . Micrometeorological convention is used here, where negative NEE values indicate ecosystem uptake of CO 2 . Sensible heat (H; W m 22 ) was determined from the covariance of the turbulent fluctuations of w and T s , and w and q (noted as primes, rf. [29]) estimated over a 30-min averaging period (noted as overbar), such that, where, r air is the air density (kg m 23 ) and C p is the specific heat of air at constant pressure (J kg 21˚C21 ). Corrections for the effect of water vapor on the speed of sound were applied [30]. Similarly, latent energy (LE; W m 22 ) was calculated from the covariance of the turbulent fluctuations of w and r v and averaged over 30-min, where, l is the heat of vaporization (J g 21 ), and M air and M w are the molecular weights of air (28.965 g mol 21 ) and water (18.01 g mol 21 ), respectively. Corrections for thermal and pressure related expansion and/or contraction, and water dilution were applied [31]. 10 Hz raw flux data were processed with EdiRe (v. 1.4.3.1184, [32]), which included despiking and air density corrections [31,33]. Fluxes (NEE, H, LE) were then corrected for mass transfer resulting from changes in density not accounted for by the IRGA [31,34], and barometric pressure data were used to correct the fluxes to standard atmospheric pressure. All measurements were filtered when systematic errors in either NEE, H or LE were indicated, such as: (1) evidence of rainfall, condensation, or bird fouling in the sampling path of the IRGA or sonic anemometer, (2) incomplete half-hour datasets during system calibration or maintenance, (3) poor coupling of the canopy with the external atmospheric conditions, as defined by the friction velocity, u*, using a threshold ,0.15 m s 21 [35,36], and (4) 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 for implausible H (,2100 or .800 Wm 22 ), LE (,2100 or .800 Wm 22 ), and NEE (i.e., NEE ,230 or .30 mmol m 22 s 21 ) values, stationarity criteria, and integral turbulent statistics [37,38]. At TS, 38% and 77% of the day and nighttime data were removed, respectively. At SRS, 34% of daytime data and 70% of nighttime data were removed. Missing H and LE were then gap-filled using the linear relationship between H or LE and R n on a monthly basis. When R 2 values were less than 70%, annual relationships between R n and H or LE were used to gap fill data in that month.
Missing half hourly NEE data were gap-filled using separate functions for day and night. When PAR was $10 W m 22 , NEE data was gap-filled using a Michaelis-Menton approach (NEE day ; Eq. 6), and when PAR was ,10 W m -2 , NEE data was gap-filled using an Arrhenius approach (NEE night ; Eq. 7): where, a is the apparent quantum efficiency ( CO 2 Lw ), w is PAR, R eco is ecosystem respiration (mmol CO 2 m 22 s 21 ), and P max is the maximum ecosystem CO 2 uptake rate (mmol CO 2 m 22 s 21 ).

NEE night~Reco~R0 exp b{T air ð7Þ
where, R 0 is the base respiration rate when air temperature is 0˚C and b is an empirical coefficient. In equation 6, R eco is an estimated model parameter, whereas R eco measurements are the dependent variable in equation 7. A bootstrap method was used for error estimation of gap-filled values of NEE. Synthetic datasets (1000) of size n (with replacement) from the original dataset of size n for each estimated gap-filling model (Eq.6 and Eq.7) on a monthly or annual basis where appropriate were used [11]. Distributions of each model parameter were constructed, which were then checked to ensure that the model parameters derived from the original data were contained within a 95% confidence region. Following gap filling, GEE was calculated from half hourly NEE and R eco data (Eq. 8).

Defining Seasons
Although the majority of rain in the Everglades region falls in the wet season, it is difficult to identify the exact onset of the wet season. Previous studies define season based on the calendar year [39] or water levels [11,12]; however, these approaches either do not capture interannual variations or are heavily influenced by water management activities performed by the South Florida Water Management District. To determine the date of the shift in seasons we examined the seasonal pattern of Bowen ratios over time, where, the subscript t denotes the t th daily value in the time series. Similar to methods used by Nuttle [40] to define hydroperiods, a harmonic function (sine function) was fit to the b time series to identify inflection points that indicate changes in the seasonal trend of the ratio of energy dissipation as H and LE (Fig. 2). A sine function was fit to the b time series at each site annually (Jan 1 to Dec 31), and the inflection point along the positively sloped portion of each sine function was used to identify the change from dry to wet season (Fig. 2). The wet season was defined by fitting the sine function to the same set of site-specific series offset by 2182 days (,6-months) and identifying the inflection point along the negatively sloped portion of each sine function. The sine function was offset by 2182 days so that the shifts in season would not occur near the end of the time series. Previous studies show marked seasonal shifts in energy dissipation in short and long hydroperiod marsh ecosystems (S. Malone, unpublished data) and therefore this method of seasonal classification should adequately capture the seasonality in both water and energy availability.
As an indication of season intensity, the seasonal mean Palmer Drought Severity Index was used (PDSI; Fig. 3a) [41]. PDSI compares weather conditions to historical weather data, taking into account temperature, rainfall, and the local available water content of the soil. PDSI uses 0 to identify normal conditions, negative numbers (21 to 26) to indicate dryer than average conditions, and positive numbers to reflect excess rain (Fig. 3). PDSI data were retrieved from the National Climatic Data Center (http://www.ncdc.noaa.gov/temp-and-precip/ drought/historical-palmers.php).

Long-term weather data
Long-term weather data were obtained from the nearest weather station, NCDC Royal Palm Ranger Station (25˚239N/80˚369W), where NOAA surface meteorological data was available from 1964 to 2013. The Oceanic Niño Index (ONI) was used to define extreme ENSO phases and was retrieved from the National Oceanic and Atmosphere Administrations Earth Physical System Research Laboratory (http://www.esrl.noaa.gov/psd/data/climateindices). The ONI is the running 3month mean sea surface temperature (SST) anomaly from a 30-year mean SST for the Niño 3.4 region (i.e., 5˚N-5˚S, 120˚-170˚W). Cold (warm) phases are defined as 5 consecutive months at or below (above) the -0.5˚(+0.5˚) anomaly.

Data Analysis
Long-term patterns in monthly weather data An intervention time series approach was used [42], utilizing autoregressive integrated moving average (ARIMA) models and the SAS procedure PROC ARIMA (version 9.3, SAS Institute Inc., Cary, North Carolina) to model 3 variables describing weather for the site, (e.g., monthly precipitation, average maximum daily temperature, and average minimum daily temperature) as a function of ENSO phase. ARIMA models incorporate 3 types of processes: autoregressive (AR) of order p, moving average (MA) of order q, and if necessary, differencing of degree d [42]. ARIMA models fit to time series data use AR and MA terms to describe the serial dependence, and use other time series data from independent variables to explain the dependence on outside factors [42]. The advantage of using ARIMA models is that the internal structure of the data (e.g., autocorrelation, seasonality) is explicitly accounted for by incorporating past values. For example, a first-order autoregressive moving average, or ARIMA(1,1) model, predicts the current time period value (Y t ) using its one period previous value (Y t-1 ) and its associated error (e t-1 ): where: m is the mean of the series, a 1 is the first-order AR coefficient, h is the firstorder MA coefficient, e t is the current period error, and e t-1 is its one period previous error. ''Memory'' can be added to the model by adding lags and associated MA and/or AR components, and potential covariates can also be added as predictor variables, which may also have lagged components.
To facilitate the inclusion of independent categorical variables for ENSO, El Niño and La Niña phases were coded as indicator variables, where a value 0 or 1 specified the absence or presence, respectively, of a categorical effect. We then determined if there were teleconnection lags between ENSO phase and precipitation and temperature.
In developing time series models, first, all data series were tested for stationarity via the augmented Dickey-Fuller unit root test [43]. The ARIMA models were then fit to time series data (monthly precipitation, and minimum and maximum temperature) using an iterative Box-Jenkins approach, where: (1) autocorrelation and partial autocorrelation analysis were used to determine if AR and/or MA terms were necessary for the given time series, (2) model coefficients were calculated using maximum likelihood techniques and, (3) autocorrelation plots of model residuals were examined to further determine the structure of the model [42].
Because of the presence of autocorrelation in the explanatory series, input series were pre-whitened [42]. ARIMA models were then fit to the dependent variables using the pre-whitened explanatory series as predictor variables. Cross-correlation coefficient plots between the explanatory series and dependent variables were used to identify direct and inverse relationships at various lags or time shifts, and autocorrelation plots of the residuals verified that the residual series had characteristics of random error, or white noise [42]. Model selection was based on minimum Akaike's information criterion (AIC) and models were acceptable when residual white noise was minimized [44]. A backwards selection method was used, removing the least significant parameter one at a time until all regression terms in the final model were significant at the a50.05 level and the lowest AIC was achieved. ARIMA assumptions of normality and independence of residuals [42] were verified by examining residual plots.

Seasonal light and temperature response
To examine changes in the response of NEE and R eco to light and temperature, respectively, non-linear equations (6) and (7) were fit. Parameters for these models were fit by ENSO phase, site and season via the SAS procedure PROC NLIN (version 9.3, SAS Institute Inc., Cary, North Carolina). Parameter estimates were then compared to identify differences and similarities. As a result of the high degree of autocorrelation inherent in NEE time series from half-hourly data, the standard errors of parameter estimates from these models are artificially small, and statistical tests are not valid. Therefore, this analysis is presented in a descriptive context.

Daily CO 2 and water dynamics
An intervention time series approach was used to identify and model the relationship between CO 2 dynamics (NEE, GEE, and R eco ) and a set of explanatory variables over a 5-year time series of daily data (2009 to 2013). These variables included: water level, season, ENSO phase (El Niño, La Niña and neutral), daily precipitation, drought condition, and average air temperature. The combined effect of ENSO phase and season (e.g. El Niño 6 wet season and El Niño 6 dry season) on CO 2 fluxes were included in time series models as predictors with indicator variables. Indicator variables were also developed to identify sections of each season that directly followed an ENSO phase (post-La Niña and post-El Niño). Beckage et al. [3] found the effect of extreme ENSO phases during seasonal transitions, and post-La Niña and post-El Niño phases capture the transition periods. To explore the effect of precipitation on CO 2 exchange rates, indicator variables were used to identify the day of a precipitation event (Rain Day), the day after precipitation, and the quantity of precipitation (Rain; mm). The indicator for the day after a precipitation event identified the first rain free day following a day with precipitation. Finally, drought conditions were defined as those days where PDSI,-2 and verified the drought extent with the National Oceanic and Atmospheric Administration's Drought Monitor ( Fig. 3b; [45]). Drought Monitor data was obtained from the National Drought Monitor Center (http:// droughtmonitor.unl.edu/DataArchive.aspx).
In addition, previous studies identified water level as one of the most important drivers for CO 2 exchange rates in the Everglades freshwater ecosystems [11,16]. Thus, a water index equal to the difference between each half-hourly water level and its site-specific annual seasonal mean water level was computed. Using the water index as a dependent variable in an additional analysis, time series models were estimated to answer questions about the relationship between intraseasonal fluctuations in water levels, precipitation, and ENSO phase.
As in the models of monthly weather data, all daily time series were tested for stationarity and non-stationary series were made stationary by differencing [46]. ARIMA models were then fit to time series data (NEE, GEE, R eco , and the Water Index) using pre-whitened explanatory series as predictor variables. Crosscorrelation coefficient plots identified relationships at various lags, and autocorrelation plots were used to verify that the residuals had characteristics of random error. Model selection was based on minimum AIC, removing the least significant parameters. For each dependent variable, a single model form was selected with common predictor variables to aid site comparisons. Non-significant parameters remained in the model only if the parameter was significant at one site and it did not affect the final model of the other site. Multicollinearity between explanatory variables was also explored to ensure models did not contain input series that were highly correlated.

Long-term weather patterns
The long term ONI data ranged from 2.5 to 22.05 and indicated that the neutral phase occurred over about half (46%) of the period 1964-2013, while the El Niño and La Niña phases occurred 26 and 27% of this long-term study period, respectively. Neutral phases ranged in length from ,1 month to.12 months and just 6 years contained no neutral phase (1969, 1971, 1975, 1987, 1999, and 2000). There were 14 different El Niño and 14 La Niña events, and the average ONI index was ,1.1, 21, and 20.02 for all El Niño, La Niña and neutral phases, respectively. Although it is common for an El Niño event to be separated from a La Niña event by a short neutral phase (14 occurrences), a neutral phase occurring between consecutive El Niño (1 occurrence) or La Niña phases (4 occurrences) were less frequent. The majority (,52%) of wet season months (May to October) were associated with a neutral phase while just ,25% and ,23% of wet season months occurred during El Niño and La Niña phases, respectively. Dry season months (September to April) were associated with the neutral phase ,41% of the time, while El Niño and La Niña phases occurred 27 and 32% of time. Time series analysis of long-term monthly precipitation and minimum and maximum daily temperatures versus ENSO phase showed that rain increased the month after the start of El Niño phases (i.e., at a lag of 1 month; p50.1043) and declined the month following the start of La Niña (p50.7719), though not significantly (Table 1; S1 Fig.). Monthly average maximum daily temperatures were lower the month following the start of La Niña (p50.0001), and during El Niño phases (p,0.0001) and at a lag of 1 month (p50.001), compared to neutral phases. Average minimum daily temperatures were lower the month following the start of El Niño (p50.019) and La Niña (p50.0218; Table 1) (Fig. 3a). During La Niña phases, seasonal mean PDSI ranged from ,21 to 24 (Fig. 3a). In the wet season of 2009 an El Niño phase began soon after a La Niña phase. The El Niño phase extended into the dry season of 2010 where it co-occurred with wetter than average conditions (PDSI.1). During the neutral phase in 2013, mean wet season PDSI was 1.7 and 20.6 in the dry season.
Although TS and SRS had similar weather, harmonic analysis of b (Fig. 2) showed that both hydroperiods and season length differed annually and between sites (Fig. 3b). The onset of the wet season at TS lagged SRS by approximately one month on average, with wet season length varying between 179 to 208 days at TS and 159 to 242 days at SRS. Wet season length was positively correlated with cumulative precipitation from January to March (p50.1495; Fig. 4a). During abnormally dry years (PDSI,22), the wet season was shortened by about 15 days at TS and 34 days at SRS compared to all other years. In 2009 and 2011, south Florida experienced severe drought conditions (water levels below the soil surface) resulting in 65 and 34 dry days at TS and SRS, respectively (Figs. 3 and 4). Drought years were not characterized by lower annual precipitation but by lower rainfall and fewer rain events the 3 months prior to the start of the wet season, which generated a shorter season (Fig. 4a). In 2010, 2012 and 2013, total rainfall

Annual and Seasonal patterns in CO 2 Fluxes
Annual net CO 2 exchange rates at both freshwater marsh ecosystems co-varied with ENSO phase (2009-2013), which corresponded to changes in season intensity. The ratio of GEE to R eco suggests that CO 2 exchange rates at both sites were most similar in 2010 (El Niño; S2 Fig.) when hydraulic conditions were comparable between TS and SRS. The lowest rates of NEE (greatest CO 2 uptake) occurred during the La Niña phase at TS followed by the neutral phase. At SRS the greatest CO 2 uptake rates were during the El Niño and neutral phase compared to seasons associated with a La Niña phase ( Table 2). The annual ratio of GEE to R eco and annual rates of NEE both show that sites were most similar when water levels were higher and during the neutral phase (S2c Fig.). Although TS ranged from a small CO 2 sink to a small source on an annual basis over the 5 years, TS was often a source for CO 2 during the wet season and a sink during the dry season (Table 2; Figs. 5 and 6), except in 2011 and 2013 when TS was a sink in both seasons. Changes in GEE relative to R eco resulted in seasonal shifts in NEE, though there was no consistent pattern in dry season versus wet season response in GEE or R eco . The ratio of GEE to R eco showed that in 2011 both seasons were comparable at TS and they differed the most in 2012 when a portion of the dry season was during a La Niña phase and the wet season was in a neutral phase (Table 2; S2a Fig.). CO 2 uptake rates were generally higher in the dry season and during the exceptionally dry La Niña years at TS, which corresponded to  (Table 2; Fig. 5c). The mean annual dry season length at TS was 172 days for the 5 years. During years with La Niña phases, CO 2 uptake was higher, dry seasons were 10 days longer on average, and TS was a greater sink for CO 2 ( Table 2). As a result of an extended drought in 2011 that occurred with 2 La Niña phases, TS was a sink for CO 2 in both wet and dry seasons ( Table 2).
Like patterns observed at TS, annual variation in NEE corresponded to changes in ENSO phase at SRS. Although SRS ranged from CO 2 neutral to a small source of CO 2 to the atmosphere annually and seasonally, CO 2 release rates increased during seasons with La Niña phases. SRS was a small source of CO 2 in most seasons, except in 2010 (El Niño) and the dry season of 2013 (neutral phase) (S2b Fig.).
Ecosystem respiration was the primary control on annual ecosystem carbon balance ( Table 2) and dry season mean R eco was often higher than wet season R eco , increasing CO 2 release at SRS when the dry season was extended (S2b Fig.). Similar to TS, average daily CO 2 uptake rates were higher in the dry season at SRS ( Table 2). Even so, dry season mean daily R eco rates were also higher than in the wet season (Table 2). Because R eco increased relative to CO 2 uptake in the dry season ( Table 2; S2b Fig.), longer dry seasons were associated with greater CO 2 source status at SRS. During the exceptionally wet year that corresponded to an El Niño phase (2010), GEE surpassed R eco and the site became a larger sink for CO 2 ( Table 2). Although site differences were apparent in patterns of dry season CO 2 exchange rates, there was no clear pattern in wet season CO 2 uptake rates at TS and SRS (Table 2; S2Fig. a and b; S3 Fig. b and d). At TS wet season NEE during the neutral phase was most similar to that of the La Niña phase at TS (Table 2; S3 Fig. a and b), while at SRS all wet seasons were similar except the 2011 wet season which co-occurred with the strongest La Niña phase ( Table 2; S3c Fig.). Plots of dry season length versus dry season cumulative NEE ( Fig. 4b and 4c) revealed differences by site. Although there were very few observations available since the study period included just 5 dry seasons, NEE at TS exhibited a negative linear relationship with dry season length (p50.0393; Fig. 4b), while NEE at SRS showed a positive linear relationship (p50.0206; Fig. 4c). These results demonstrate that the seasonal response in NEE rates differed between sites, and suggest that dry season length (and changes in dry season length) may control the CO 2 source and sink status in the future. Although, the data set is not yet large enough to confirm the pattern between dry season length and dry season cumulative NEE as a characteristic of each site, this relationship is an important indication of how the sites respond differently to hydroperiods.

Seasonal light and temperature response
In addition to site differences, ENSO phase and season altered photosynthetic capacity and ecosystem respiration (Table 3; Fig. 6). At both sites, photosynthetic capacity (P max ) was greatest during La Niña and directly following La Niña phases compared to rates during neutral and El Niño phases, although at higher PAR values (.1500 mmol m 22 s 21 ) curves for La Niña and the neutral phase converged (Fig. 6a and 6b). Seasonal differences in light and temperature response curves were greater, and photosynthetic capacity and dark respiration (R eco ) were consistently higher at TS, for all ENSO phases as compared to SRS (Table 3; Figure 6a and 6c). At both sites, the effect of El Niño and La Niña phase was also greater during the dry season. At SRS, photosynthetic capacity was similar for all ENSO phases during the dry season though the effect of El Niño and La Niña phases increased during the wet season (Table 3; Fig. 6b and 6d;). Overall, there was a small seasonal difference in photosynthetic capacity, which was higher on average in the wet season than in the dry season at SRS. At both sites, differences in photosynthetic capacity by ENSO phase were greatest at higher PAR values (.1000 mmol m 22 s 21 ). The relationship between temperature and R eco differed between ENSO phases, and R eco showed distinct seasonal patterns in temperature sensitivity (Table 3; Fig. 6e to 6h) at TS and SRS. During the wet season at both sites, R eco was less sensitive to temperature changes. At TS temperature effects associated with ENSO phases were greater at lower temperatures (Table 3; Fig. 6e and 6f). At higher temperatures R eco was more sensitive to changes in temperature during all phases at both sites, a response that was enhanced during the dry season (Table 3; Fig. 6g and 6h). At SRS, the differences among ENSO phases were small except at high temperatures (.24˚C) during the wet season, while at TS the differences among ENSO phases were consistently large at lower temperatures and converged at high temperatures ( Fig. 6e to 6h). Similar to patterns observed in light response curves, respiration rates were higher at TS than at SRS, and temperature patterns associated with R eco also showed greater release of CO 2 at both sites in the dry season versus the wet season (Fig. 6e to 6h).
The effect of ENSO phase, precipitation, and season on daily CO 2 exchange rates and water level After pre-whitening, some small (,0.05) but statistically significant autocorrelation remained in pre-whitened series; however, this sensitivity resulted from the large number of observations available and was judged to be biologically insignificant [47]. Differencing was required for water level and PDSI time series due to the lack of stationarity at both sites. Non-stationarity indicates a lack of stability in the mean of these variables over time, further suggesting that there were significant changes in hydroperiods at both sites. By including differenced variables in models we evaluated how changes in water level and PDSI influenced NEE, R eco , and GEE. In response to evidence of 1-month lagged teleconnections for ENSO phase effects on precipitation and average daily maximum and minimum temperatures (Table 1), lagged El Niño and La Niña phase indicators were included in time series models of CO 2 exchange rates. Models for NEE (S4a and S5a Figs.), had a significant lagged 1-day MA [MA(1)] component, as well as significant AR components at a lag of 1-, and 2days (p,0.0001; Table 4). At TS, Dwater level (p50.01) and the quantity of rain (mm day 21 ; p,0.0001) had significant positive relationships with NEE, showing that as the change in water level and the quantity of daily rainfall increased, net CO 2 uptake decreased (higher NEE; Table 4). Post-La Niña phases in the dry season were associated with significantly lower NEE (higher net CO 2 uptake) compared to neutral and El Niño phases at TS (p50.0039); however, post-La Niña phases in the dry season at SRS were associated with greater NEE (lower net CO 2 uptake; p50.0394). Moreover, there was a significant increase in NEE at TS (higher net CO 2 uptake) the day after rain. At TS where hydroperiods were shorter, the effect of rain and post-La Niña phase during the dry season were significantly stronger than at SRS ( Table 4). The quantity of rain was the strongest driver of NEE at both sites.
Models for R eco (S4b and S5b Figs.), at both sites had a significant 1-day MA [MA(1)], as well as significant AR components at a lag of 1-and 2-days (p,0.001; Table 5). The day after rain (p50.0047), and post-La Niña phases in the dry season (p50.0291) reduced daily ecosystem respiration rates at TS (Table 5), while the quantity of rain (p50.0596), and post-El Niño phases in the dry season were associated with an increase in R eco . At SRS, post-La Niña phases in the dry season (p50.0002) were significantly linked to increased R eco (Table 5), and the effect of post-La Niña phases in the dry season was larger at SRS where it was associated with increased R eco . At SRS, post-La Niña phases had the greatest impact on R eco , while at TS R eco had the strongest association with the day following rain events.
Similar to NEE, models for GEE had a significant 1-day lagged MA [MA(1)] and AR [AR(1)] component at both sites (p,0.0001; Table 6; S4c and S5c Figs.). The day of rain (p,0.0001), day after rain (p50.0015), quantity of rain (p,0.0001), and Dwater level (p,0.0001) had significant positive relationships with GEE at TS (Table 6). At SRS, days with precipitation (p,0.0001) and Dwater level (p50.005) had positive relationships with GEE (Table 6). Rain had a stronger effect at TS than at SRS, and rain had the strongest effect on GEE at SRS, while Dwater level had the strongest effect on GEE at TS (Table 6).
Models for the water index had a significant 1-day lagged MA [MA(1)] and AR [AR (1)] component at both sites (p,0.0001; Table 7; S4d and S5d Figs.). At TS, La Niña phases during the wet season (p,0.0001) and post-La Niña phases (wet season; p,0.0001) were associated with lower than average water levels compared to neutral and El Niño phases (Table 7). Water levels at TS were higher than the seasonal average the day of rain (0.0545), the day after rain (p50.162), during dry season EL Niño phases (p,0.0001), and throughout post-La Niña phases that occurred at the end of the dry season (p,0.0001). At SRS, the day of rain (p50.0349), the day after rain (p50.0066) and post-La Niña phase at the end of the dry season (p,0.0001) were associated with higher than average water levels ( Table 7). La Niña and post-La Niña phases (wet season; p,0.0001) were  (1) is the estimated moving average term at a 1-period lag (1 day), and AR (1) and AR(2) are the estimated autoregressive terms at a 1-and 2-period lags (1 and 2 days). Lagged values of independent variables are denoted similarly. Asterisks denote significant differences between sites. Day after rain is an indicator for the first rain free day, Rain is the quantity of precipitation (mm), DWater Level is the change in water level from one day to the next, Post-La Niñ a (Wet Season) is an indicator for the time directly following a La Niña phase in the wet season, and Post-La Niñ a (Dry Season) is an indicator for the time directly following a La Niña phase in the dry season.
doi:10.1371/journal.pone.0115058.t004  associated with lower than average water levels. The effect of La Niña (wet season) and El Niño (dry season) phases were significantly greater at TS than at SRS, and La Niña and post-La Niña phases were the strongest predictors of the water index (Table 7).

Discussion
The goal of this research was to understand the relationship between ENSO phases and CO 2 exchange rates (NEE, R eco and GEE) in Everglades freshwater marsh ecosystems. The relationships between ENSO extremes, precipitation and hydrology in the Everglades region suggest that El Niño and La Niña phases could be important for C dynamics [1,3,48,49,50,51,52,53,54,55]. The results presented here show that climate teleconnections have significant controls on Everglades CO 2 dynamics, demonstrating that in addition to climate change and water management, ENSO is an additional source of variation in C cycling in Everglades freshwater ecosystems.

Annual fluctuations in season length and intensity
Interannual variability in precipitation can be large in the Everglades region [21,56], which was the case throughout the study period. Reduced precipitation prior to the onset of the wet season resulted in shorter wet seasons and lower water levels, though annual precipitation was unchanged. Cumulative precipitation in January through March was related to wet season length (Fig. 4a).
Although more data is needed to validate this relationship, previous research suggests dry season rainfall from October to April largely determines season intensity [3,56]. Higher precipitation from January to March was also associated with longer and wetter than average wet seasons. Interannual variation in the onset and length of seasons can have a significant effect on the magnitude of ecosystem primary production [57], and changes in season intensity can either suppress or enhance production [58]. Results suggest that in Everglades ecosystems dry season length has a strong relationship with annual NEE (S6 Fig.).
Knowing that precipitation patterns were driving the variations observed in season length and intensity, we examined the co-occurrence of ENSO phases, previously found to alter season intensity in the Everglades region [3,55], to PDSI defined season intensity. Abnormal precipitation and water level patterns over the study period coincided with El Niño and La Niña phases. Changes in ENSO phase have been associated with shifts in the position of the midlatitude jet, which is important for patterns in frontal precipitation. In the Everglades where frontal precipitation is the main source of dry season rainfall [20], changes in ENSO phase can have a significant effect on hydroperiods. During El Niño (La Niña) phases, the equatorial (poleward) displacement of the midlatitude jet increases (decreases) frontal precipitation in the southeastern United States [1,48,49]. Studies have shown that in Florida, El Niño is positively correlated with winter (dry season) precipitation, explaining up to 34% of dry season precipitation variability [4]. Here, El Niño and La Niña phases have been shown to reduce seasonal differences in rainfall without altering annual precipitation inputs [55].
Results found by Beckage et al. [3] and our study show that El Niño (La Niña) phases were correlated with increased (decreased) rainfall and water levels in the Everglades region. Time series analysis of the water index at TS and SRS reinforced the relationships previously found between precipitation patterns and ENSO phase. At both TS and SRS, La Niña and post-La Niña phases were associated with lower than average water levels and at TS water levels were higher than the seasonal average during the El Niño phase (dry season). The effect of La Niña (wet season) and El Niño (dry season) phases were significantly greater at TS than at SRS, and La Niña and post-La Niña phases were the strongest predictors of the water index at both sites. These results support the previously observed patterns in ENSO phases [3], precipitation, and hydrology and show that El Niño and La Niña phases induced fluctuations in season intensity (+ and -, respectively) compared to neutral phases, which affect season length, and have important implications for annual CO 2 exchange rates.

Seasonal patterns in CO 2 exchange rates
Hydroperiods have shaped soil conditions and species composition at each site in ways that have led to different seasonal patterns in CO 2 exchange rates. Hydroperiods alter ecosystem production by interfering with exposed leaf area [11,12,16], triggering senescence (S. Oberbauer, unpublished data), and allowing CO 2 fixation within the water column [16]. Seasonal changes in photosynthetic capacity (Fig. 6) and the relationship between respiration rates and temperature support the patterns previously found in Everglade freshwater marsh studies [11,12]. Knowing that season intensity changed with El Niño (+) and La Niña (-) phases (Fig. 3a), we expected and saw a magnification of the site-specific seasonal response in CO 2 exchange rates during and directly following El Niño and La Niña phases. The lag in El Niño and La Niña phase effect on CO 2 exchange rates is the result of the effect of extreme ENSO phases on water levels during transition periods. El Niño (La Niña) phases increase (decrease) surface water levels during seasonal transitions and, major drainages often contain no water during transitions in a La Niña phase [3]. Results suggest that Everglades freshwater marsh ecosystems are more similar during El Niño and neutral phases when water levels are higher. This is further reinforced by time series analysis, which detected significant changes in CO 2 exchange rates between ENSO and post-ENSO phases at TS and SRS.
The sites differed in their response to ENSO phases. At TS, CO 2 exchange rates were sensitive to El Niño and La Niña phase in both the wet and dry seasons, while at SRS the effect of El Niño and La Niña phase on the wet season was not as strong as the effect in the dry season. La Niña phases resulted in higher photosynthetic capacity and greater seasonal net carbon uptake rates in both seasons at TS compared to neutral and El Niño phases. While La Niña increased photosynthetic capacity at SRS, the increase in R eco reduced net exchange rates compared to neutral and El Niño phases. As a result of changes in GEE relative to R eco , El Niño led to greater net CO 2 uptake at SRS during the dry season compared to La Niña and neutral phases. As water levels and precipitation increased, NEE and GEE also increased at both sites (less CO 2 and net CO 2 uptake). At TS where hydroperiods were shorter, the effect of rain and the post-La Niña phase during the dry season were significantly stronger than at SRS (Table 4). Rain was one of the strongest drivers of CO 2 dynamics and these results indicate that there is a significant relationship between El Niño and La Niña phase and season intensity, either creating conditions wetter (+) or dryer (-) than normal, which magnifies CO 2 exchange rates at TS and SRS.

Annual patterns in CO 2 exchange rates
The results presented here show that the length and intensity of the wet and dry season varied annually with climate patterns. Considering the site-specific response to season, these results support our hypothesis that variations in season length would explain interannual fluctuations in NEE (S6 Fig.), R eco , and GEE. At TS where mean GEE surpassed R eco during the dry season, an increase in dry season length and intensity amplified the site's net CO 2 uptake rates. Annually, SRS was usually a small source of CO 2 , although when wet season conditions intensified during El Niño, net CO 2 uptake increased. The effect of ENSO phase also differed by site, showing that longer hydroperiods mute the effect of climate fluctuations on CO 2 exchange rates, and as water levels decline the system becomes more vulnerable to climate. The ecosystem's sensitivity to climate fluctuations has important implications for water management and climate change [59]. The uncertainty of climate change makes it important to understand how ecosystems respond to climate events and how these responses aggregate to form trends in net CO 2 exchange rates.

Effect of climate change and water management on net CO 2 exchange rates
In sub-tropical ecosystems, phenology is less sensitive to temperature and photoperiod, and more tuned to seasonal shifts in precipitation [60,61,62]. Such shifts are expected to occur in concert with rising global temperatures, but both the direction and magnitude of change vary regionally [58,63]. As climate change has the potential to alter hydrologic regimes, we can expect to see greater variations in CO 2 exchange rates. Shifts in water management and land use change (e.g., conversion to agriculture and urban development) could also significantly alter both hydrology and CO 2 dynamics in the Everglades region, making it important to develop a baseline understanding of how hydroperiod drives changes in CO 2 dynamics and how climate alters hydrology. With water managers striving to adjust hydroperiods closer to natural values, in the future we might expect water levels in TS to increase, offsetting changes in climate by maintaining current patterns in hydrology. Alternatively, we might anticipate higher water levels to increase hydroperiods, making the system less sensitive to climate change altogether. With longer hydroperiods, SRS will likely remain less sensitive to changes in climate and land management. As a result, the two sites will likely behave more similarly in the future as SRS remains neutral or a very small source of CO 2 to the atmosphere and TS becomes more neutral.
Patterns observed in ENSO phases and the co-occurrence of extreme wet and dry seasons suggest changes in climate patterns can significantly alter ecosystem function. Equatorial Pacific SST during the past half century show a clear warming trend that is consistent with global warming [64], and El Niño and La Niña phases are expected to continue increasing in severity and frequency [3,17,58]. Moreover, as a result of warming SST, ENSO amplitude may become even stronger, intensifying feedbacks relevant to ENSO phases [64]. If the frequency and intensity of strong climatic disturbances increases beyond historical averages, altered disturbance regimes have the capacity to significantly modify ecosystem processes [65]. ENSO phases have been linked to climate anomalies [53] and CO 2 dynamics [66] on a global scale [53], making it crucial to analyze the importance of ENSO extremes and other cyclic climatic phenomenon on the variability of terrestrial carbon cycling.

Study Limitations
This research shows that the length and intensity of the wet and dry season vary annually with ENSO phase in the subtropical Everglades region where dry season precipitation is dependent on frontal systems. Changes in season length and intensity are also correlated with CO 2 exchange rates showing that extreme ENSO phases magnify the site-specific seasonal response in freshwater marsh ecosystems. Although a longer time series is required to verify that this relationship is persistent, this research provides initial insights into an important driver of seasonal and interannual variation in CO 2 exchange rates in Everglades ecosystems.