Earlier snowmelt may lead to late season declines in plant productivity and carbon sequestration in Arctic tundra ecosystems

Arctic warming is affecting snow cover and soil hydrology, with consequences for carbon sequestration in tundra ecosystems. The scarcity of observations in the Arctic has limited our understanding of the impact of covarying environmental drivers on the carbon balance of tundra ecosystems. In this study, we address some of these uncertainties through a novel record of 119 site-years of summer data from eddy covariance towers representing dominant tundra vegetation types located on continuous permafrost in the Arctic. Here we found that earlier snowmelt was associated with more tundra net CO2 sequestration and higher gross primary productivity (GPP) only in June and July, but with lower net carbon sequestration and lower GPP in August. Although higher evapotranspiration (ET) can result in soil drying with the progression of the summer, we did not find significantly lower soil moisture with earlier snowmelt, nor evidence that water stress affected GPP in the late growing season. Our results suggest that the expected increased CO2 sequestration arising from Arctic warming and the associated increase in growing season length may not materialize if tundra ecosystems are not able to continue sequestering CO2 later in the season.

www.nature.com/scientificreports/ mental controls during different periods of the summer corresponding to various stages in seasonal phenology (early season: June, peak season: July, and late season: August). We used a partial correlation analysis to identify if the timing of the snowmelt associates with anomalies of ET, soil moisture, NEE, GPP, ER, atmospheric vapor pressure deficit (VPD), or the Bowen ratio (the ratio between Sensible Heat (H) and Latent Heat (LE)) while statistically controlling for the main meteorological forcing such as air temperature and solar radiation (Methods). Identifying the correlation between ET (and the Bowen ratio) and snowmelt timing is a way to assess water limitation to ecosystems (in addition to testing their response to soil moisture changes), as H, and therefore, the Bowen ratio, are expected to increase with surface drying 31,32 . To identify the association between snowmelt timing, the main environmental variables (i.e., air temperature and solar radiation), and NEE, GPP, ER, and ET over time, we performed a maximum covariance analysis (MCA) on the monthly median standardized anomalies from 2004 to 2019 (a time period when data for most of the sites were available). MCA allowed us to find patterns in two space-time datasets that are highly correlated using a cross-covariance matrix 26 . We retained sites as the unit of variation (i.e., by estimating the standardized anomalies by site for each of the indicated variables, see "Methods"), hence the results of the MCA integrated the site level relationships between each of the variables over time). The goal of this analysis was to identify the most important environmental drivers associated with NEE, GPP, and ER across all the sites over time. MCA is particularly appropriate for this study as it can handle data with gaps and unequal lengths in the datasets. We also tested the relative importance of the abovementioned environmental drivers on the monthly median GPP, ER, and NEE using a linear mixed effect model, including site as a random effect to account for the site-to-site variability. The MCA and the mixed model analyses were conducted to test the relative importance of snowmelt and other variables at different times of the season. Finally, to evaluate the water balance at different times of the season, we estimated the difference between Potential Evapotranspiration (PET) and the actual ET, and the difference between precipitation (PPT) and ET for each of the sites, years, and months (e.g. June, July, and August). This study did not attempt to describe the long-term temporal changes in the anomalies of snowmelt and carbon fluxes, given the short data record available for some of the sites (i.e. less than 10 years, Table S1), but instead focused on understanding the association between environmental variables and the carbon balance at different times of the season. More details of these analyses are included in the Methods.
Influence of snowmelt timing on NEE, GPP, ER, and hydrological status of tundra ecosystems. Once statistically controlling for solar radiation and air temperature (in the partial correlation analysis, see "Methods"), we observed a significant positive relationship between the snowmelt timing anomalies and NEE anomalies (i.e. earlier snowmelt was associated with a higher net CO 2 sequestration) in June and July, but a negative correlation in August (Fig. 2a, Table 1). A significant relationship was also found between snowmelt date anomalies and GPP anomalies, with more positive GPP anomalies (i.e. higher plant productivity) with earlier snowmelt in June and July, and more negative GPP anomalies with earlier snowmelt in August (Fig. 2b, Table 1). Earlier snowmelt was associated with significantly higher ER in both June and July, but there was no significant relationship in August (Fig. 2c, Table 1), suggesting that the late-season correlation between NEE and snowmelt timing was mostly driven by the lower GPP and with earlier snowmelt in August. The MCA analysis showed that  www.nature.com/scientificreports/ the anomalies in snowmelt timing had the highest squared covariance fraction (SCF) with the monthly median anomalies of GPP, NEE, and ER in June and July, and the lowest in August over the 2004-2019 period (Fig. 3). A similar result was observed in the linear mixed effect model, which showed a significant relationship between snowmelt date and GPP, and NEE, in all summer months, higher R 2 m between the snowmelt date and GPP in June and July, and no significant relationship between snowmelt date and ER in August (Table S3). In late season, other environmental variables had a higher covariance with the GPP, NEE, and ER anomalies than the snowmelt timing (Fig. 3, Table S3).
Our results are consistent with the discrepancy between the observed increase in the maxNDVI over the last four decades and the time-integrated (TI) NDVI which instead has plateaued in the last two decades and even decreased over the last 10 years in several northern arctic ecosystems 33 . TI-NDVI considers the length of the growing season and phenological variations 34 and, therefore, better integrates vegetation development during the entire growing season. Moisture was shown to be important for the NDVI trends 33,35 . Given the potential water limitation to summer carbon uptake in northern ecosystems 12,23-25 , we tested if an earlier snowmelt was associated with a decrease in soil moisture, which would affect GPP and NEE. We only observed a significant correlation between soil moisture anomalies and snowmelt date anomalies in June (i.e. higher soil moisture with earlier snowmelt, Fig. S1a, Table S2), but no significant correlation in July and August (Fig. S1a, Table S2). The higher soil moisture with earlier snowmelt in June is consistent with surface inundation after snowmelt 36,37 and earlier soil thawing resulting in higher soil moisture (i.e., soil moisture is low while soils are frozen). A similar result was observed for the ET anomalies. Higher ET with earlier snowmelt in June (Fig. S1b) could be the result of surface inundation after snowmelt 32 . The standardized NEE anomalies were significantly correlated with the soil moisture anomalies in each of the summer months (Fig. S1d, Table S2). However, the relationship between Figure 2. Relationships between the indicated median monthly anomalies using partial correlation analysis accounting for solar radiation and air temperature anomalies (retaining site as the unit of variation). Given that the interaction term between "month" and snowmelt timing was significant, we included the correlation coefficients and P of the regressions for each of the indicated months separately in each panel (also included in Table 1). Negative values indicate CO 2 uptake and positive values CO 2 release into the atmosphere, and shaded areas are 95% confidence intervals. Table 1. Significance (P) and Pearson's correlation coefficient (r) of the relationships between the indicated monthly median standardized anomalies for June, July, and August retaining site as a unit of variation. Site was retained as the unit of variation by estimating the standardized anomalies by month and site for each of the indicated variables. The anomalies of the indicated variables were regressed with snow depth anomalies using a partial correlation accounting for the anomalies of solar radiation and air temperature, as shown in Fig. 2. The P and r were only included when the P < 0.1 (given that for P > 0.1 we assumed that r is not significantly different from zero), N = 284. www.nature.com/scientificreports/ the GPP (and ER anomalies) and soil moisture anomalies was only significant in June (Fig. S1e,f, Table S2) suggesting an earlier activation of the vegetation with earlier soil thaw (and the associated higher soil moisture). A higher water loss from ET in early season (Fig. S1b) could have resulted in the drying of the surface moss layer with the progression of the summer, which would have been consistent with the observed lower GPP and the lower net CO 2 sequestration with earlier snowmelt observed in August (Fig. 2a,b, Table 1). A potential moisture limitation to plant productivity might have been consistent also with the higher SCF of NEE, or GPP and VPD anomalies in August than in June and July (Fig. 3). However, no significant relationship between ET (or soil moisture) and snowmelt date anomalies was observed in July and August (Fig. S1a,b) contrary to what would be expected if drying occurred following earlier snowmelt. No significant relationship was found between VPD anomalies and snowmelt date anomalies in any of the summer months (P = 0.14 in a partial correlation considering air temperature and solar radiation anomalies). Finally, surface drying should result in an increase in the Bowen ratio anomalies with the progression of the summer, given that H increases with a decrease in water table and surface drying 32,38 . However, the Bowen ratio showed no correlation with the standardized snowmelt date anomalies in any of the summer months (Fig. S1c, Table S2), and presented similar values in all the summer months (Fig. S2a). The lack of correlation between the soil moisture, VPD, Bowen ratio, and snowmelt date anomalies suggests that an earlier snowmelt did not result in significant surface drying in the sites of this study. The median PET-ET and PPT-ET for all years and sites included in this analysis (Fig.S2 b,c) was slightly higher in August, similar to reports by others for the Russian arctic tundra 38,39 , further supporting a lack of soil moisture limitation in late season. Although these analyses do not consider runoff, which can be significant 21,26 ,  www.nature.com/scientificreports/ overall our results do not suggest that an earlier snowmelt resulted in a water stress that significantly limited plant productivity in these arctic ecosystems over continuous permafrost. The correlation between the anomalies in the August GPP and snowmelt timing is consistent with earlier senescence in northern plant species (e.g. Eriophorum vaginatum, a dominant species across these tundra types) compared to southern species growing in the same location in a common garden experiment 40 . The phenotypic variation was shown to be persisting for decades 41 , and ecotypes may be unable to extend their effective growth period or take advantage of a longer growing season 40 . Several studies across different plant functional types have shown that once plant growth is initiated after the snowmelt in northern ecosystems, it continues only for a fixed number of days until the occurrence of senescence [42][43][44] . Therefore, the lower GPP in August with earlier snowmelt might not be linked to water limitation on photosynthesis later in the season, but rather to an earlier senescence arising from the endogenous rhythms of growth and senescence, that plant functional types living in these extreme conditions have developed over decades. On a broader scale, earlier senescence with an earlier start of the growing season after snowmelt in northern ecosystems is also consistent with an earlier spring zero-crossing date and an earlier autumn zero-crossing date of the mean detrended seasonal CO 2 variations at Barrow, AK, USA (NOAA ESRL: https:// www. esrl. noaa. gov/ gmd/ ccgg/ obspa ck/) during 2013-2017 compared to 1980-1984 5 . The spring and autumn zero-crossing date is the time when the detrended seasonal CO 2 variations intersect the zero line in spring and autumn respectively and can be used as an indicator for the start and end of the net CO 2 uptake by vegetation 45,46 . On the other hand, NDVI measurements show both an earlier start of the season and a later end of season for 2008-2012 compared to 1982-1986 5 . The disagreement between the detrended seasonal atmospheric CO 2 concentration showing an earlier autumn zero-crossing date and the NDVI measurements showing a later end of the season has been explained by the increase in respiration in the fall 20 . The disagreement between atmospheric CO 2 concentration trends (showing an earlier autumn zero-crossing date), and NDVI (showing a later end of the season, 5 may also be explained by the challenges in using NDVI as a proxy for plant productivity in these arctic systems. The relationship between NDVI and CO 2 flux and plant productivity is highly variable and non-linear in arctic ecosystems 47 . While some arctic ecosystems have shown that NDVI was strongly correlated with GPP (explaining 75% of the variation in GPP 48 , other studies showed that NDVI was either not significantly correlated with GPP and NEE 49 or was only able to explain a minor fraction (maximum of 25%) of the variation in NEE and GPP in some arctic tundra ecosystems after accounting for the seasonal variation 50,51 .
In conclusion, earlier snowmelt was associated with greater net CO 2 uptake and higher GPP in early and peak seasons, but with less net CO 2 uptake and lower GPP later in the summer, in the studied arctic tundra ecosystems. We did not find evidence of a late-season water limitation to GPP with earlier snowmelt. Although several hypotheses can be forwarded to explain the link between snowmelt and late season declines in plant productivity and carbon uptake, the current literature does not provide a definitive explanation (schematic Fig. 4). Future studies should investigate the potential interaction of different processes explaining the response of the carbon dynamics in the Arctic to earlier snowmelt and reconstruct the temporal changes in the carbon balance from these systems. The link between the long-term changes in the CO 2 fluxes and NDVI in tundra ecosystems needs closer examination. Studies should investigate if higher NDVI is definitively associated with higher net CO 2 uptake. Greening of the Arctic might not necessarily translate into more net CO 2 uptake, as early and peak season carbon gains might be offset by a late-season CO 2 loss, and respiration might counterbalance the increase in plant productivity. A better understanding of the processes driving these temporal changes is a fundamental step in advancing our prediction of the response of the arctic CO 2 balance to changing climate.

Materials and methods
Site description. A total of 11 eddy covariance flux tower sites across the Arctic were used in this study, where each site had at least six summers of flux data available (SI , Table S1). Ecosystem-scale CO 2 fluxes were estimated using the eddy covariance method [52][53][54] . Details pertaining to the sites, data processing, and gap-filling are provided in the SI Appendix, Table S1. All sites are located in continuous permafrost tundra regions. Vegetation at the tower sites, the instruments used to measure fluxes, the average environmental conditions at each site, the datasets used in this study for each site, and the references describing the sites are indicated in SI Appendix, Table S1. As shown in Fig. S2, the Bowen ratio reported for this study showed similar values to what previously reported during the growing season months in the Arctic (from 3.9 in a dry heath to 1.6 in a wet fen in Greenland 31 ; 0.83 38 to 0.20-0.25 in two Siberian Arctic sites 39 ; and 0.51-1.69 in a moist-tussock tundra in Alaska 32 . To estimate the standardized anomalies in the soil moisture we selected the most consistent depths and sensors (i.e., the same sensor available for the entire time period in each site, or sensors at the most similar depths in each site and across sites when data from the same sensor was not available due to instrument failure).  Supplementary Table S1.
The R package 'Evapotranspiration' (Version 1.15 55 ) was used to estimate the daily aggregated Priestley-Taylor potential evaporation 56,57 at each of the study sites, then summed into monthly totals. We subtracted monthly total actual ET measured with eddy covariance at the respective sites to estimate the PET-ET shown in Fig. S2b. Raster files of monthly precipitation accumulation were acquired for the months of June-August from www.nature.com/scientificreports/ TerraClimate 58 over the years 1959-2019. Precipitation data was then extracted for the Eddy Covariance tower coordinates using the terra package 59 in R 60 to estimate monthly total precipitation for June, July, and August for each of the sites (Tables S1). We did not use the precipitation collected by the meteorological sensors installed in the tower sites given large gaps in the site-level dataset. The other environmental variables used in this study were collected at the tower sites. The median difference between the total precipitation and the total ET in each site was estimated to evaluate the PPT-ET at each study site during each time of the season (June, July, and August as shown in Fig. S2c). The median was used as it is less affected by outliers.

Statistical analysis.
Site-level data. For the analyses performed in this study, we separated the data into different times of the season (early season: June, peak season: July, late season: August), given that some of the environmental controls could be very different given the distinct stages of vegetation development. A partial correlation analysis was carried out to identify the correlation between the monthly median standardized anomalies (the ratio between the anomalies and the climatological standard deviation) of NEE, GPP, ER, ET, snowmelt date, and other environmental variables (most of which covary). The standardized monthly median anomalies were estimated by calculating the monthly median of the standardized weekly anomalies. The standardized weekly anomalies were estimated by subtracting the observed values (for each week, year, and site) from the respective average values (for each week, and site and for the entire period available in each site), and then standardized by dividing the anomalies by the standard deviation for each week, and site for the entire period available in each site. The standardized median monthly anomalies were used instead of weekly values, given that the median is less affected by outliers. The monthly scale was used for all the variables included in these analyses as the most appropriate temporal scale to identify the importance of the variability in soil moisture on CO 2 fluxes (given that soil moisture does not change much at the hourly and weekly scale at these tundra sites). The NEE, GPP, ER data used in these analyses were gap-filled using standard methodologies as described in the Supplementary Information. The partial correlations tested the relationships between the standardized snowmelt date anomalies and the monthly median standardized anomalies of NEE, GPP, ER, ET, VPD, and soil moisture, retaining site as the unit of variation (while controlling for solar radiation and air temperature anomalies as the main controls on carbon fluxes). We also tested if the inclusion of a given site within a linear mixed model changed the results of the correlation analysis between the anomalies. To this purpose, linear mixed effects models (nlme package in R 60 ) were used to test the significance of the correlation between the above-mentioned anomalies, by including "site" as categorical random effects to account for pseudo-replication due to the different sites measured in different years. Model performance was evaluated based on the Akaike information criterion (AIC) values, on the marginal coefficient of determination ( R 2 m similar to the explanatory power of the linear  www.nature.com/scientificreports/ models) for generalized mixed-effects models as output by the "r.squaredGLMM" function within the "MuMIn" package in R 61,62 . These analyses are included in Table 1, Table S2, and Table S3.
To maximize the dataset for each analysis we included all available time periods for the variables regressed in Fig. 2, but only selected the 2004-2019 period for the Maximum Covariance Analysis (MCA) to include a time period where most sites had data available. The MCA was performed on two fields (i.e., anomalies in NEE, GPP, and ER and anomalies in the environmental drivers, see Fig. 3); the columns of the two fields are spatial locations (each site was retained as a unit of variation of this analysis) and rows are temporal measurements. Site was retained as the unit of variation by estimating the median standardized anomalies by year, week, and site for each of the indicated variables, and then by estimating a monthly median from these weekly anomalies. The first pair of singular vectors are the phase-space directions when projected that have the largest possible crosscovariance. The singular vectors describe the patterns in the anomalies that are linearly correlated. We used the time series of the first singular value decomposition (SVD) mode to visualize the parts of the datasets that vary together and reported the squared covariance fraction (SCF) of the MCA in Fig. 3. Given the limited length of the dataset we did not discuss the long-term changes in the reported anomalies. However, the MCA allowed us to evaluate the influence of snowmelt timing on the carbon balance over time at different times of the growing season, and compare its relative importance to other variables. All analyses were carried out in R 60 .