Atmospheric teleconnection influence on North American land surface phenology

Short-term forecasts of vegetation activity are currently not well constrained due largely to our lack of understanding of coupled climate-vegetation dynamics mediated by complex interactions between atmospheric teleconnection patterns. Using ecoregion-scale estimates of North American vegetation activity inferred from remote sensing (1982–2015), we examined seasonal and spatial relationships between land surface phenology and the atmospheric components of five teleconnection patterns over the tropical Pacific, north Pacific, and north Atlantic. Using a set of regression experiments, we also tested for interactions among these teleconnection patterns and assessed predictability of vegetation activity solely based on knowledge of atmospheric teleconnection indices. Autumn-to-winter composites of the Southern Oscillation Index (SOI) were strongly correlated with start of growing season timing, especially in the Pacific Northwest. The two leading modes of north Pacific variability (the Pacific-North American, PNA, and West Pacific patterns) were significantly correlated with start of growing season timing across much of southern Canada and the upper Great Lakes. Regression models based on these Pacific teleconnections were skillful predictors of spring phenology across an east-west swath of temperate and boreal North America, between 40°N–60°N. While the North Atlantic Oscillation (NAO) was not strongly correlated with start of growing season timing on its own, we found compelling evidence of widespread NAO-SOI and NAO-PNA interaction effects. These results suggest that knowledge of atmospheric conditions over the Pacific and Atlantic Oceans increases the predictability of North American spring phenology. A more robust consideration of the complexity of the atmospheric circulation system, including interactions across multiple ocean basins, is an important step towards accurate forecasts of vegetation activity.


Introduction
Photosynthetic uptake of CO 2 by terrestrial vegetation removes approximately 25% of anthropogenic CO 2 emissions from the atmosphere (Ciais et al 2013), and the magnitude (i.e. primary production) and seasonality (i.e. phenology) of this vegetation carbon uptake leave distinct fingerprints on atmospheric CO 2 concentrations (Tucker et al 1986, Keeling et al 1996, Graven et al 2013, Forkel et al 2016. However, short-term forecasting and long-term projection of vegetation activity are currently hampered by comparatively poor understanding of the drivers of vegetation change , particularly at interannual time-scales (Luo et al 2015). A better understanding of the causes of variability in land surface phenology and other ecosystem processes could improve our ability to forecast vegetation activity, which is increasingly desirable within applications ranging from agriculture to fire and land management to human health (White andNemani 2006, Morisette et al 2009). Over longer time-scales, secular trends in temperature and hydroclimate will be overlain on internal climate variability, and any increases in the variability of the climate system will exacerbate changes in mean climate state (Cubasch et al 2013, Frank et al 2015. Understanding present sources of variability and their impacts on ecosystems is therefore necessary when considering possible effects of climate change.
Much of the internal variability in the climate system arises from coupled ocean-atmosphere teleconnection patterns (Hartmann et al 2013), which affect climate on global scales, even in regions that are located far from the oceanic and atmospheric centers-of-action. Capitalizing on lags between these circulation systems and their teleconnected effects on surface climate could allow for seasonal predictions of vegetation activity in advance of the growing season. Indices describing atmospheric circulation patterns, particularly the El Niño-Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO), are already important features of many seasonal weather forecasting systems (e.g. Fereday et al 2012, Kalra et al 2013 due to their impacts on regional climates and predictive lead-times of up to several months. These teleconnection patterns are typically treated individually or independently of each other, but in reality they describe components of a complex atmospheric circulation system in which the phase of one teleconnection may modulate the climate effects of another . Understanding the predictive capability of atmospheric teleconnection patterns on ecosystem processes and capturing interactions among teleconnections could therefore potentially allow for seasonal forecasting of vegetation activity over large areas. At present, some of these circulation patterns, particularly ENSO, may be more variable than in the paleoclimate past (Cobb et al 2013, Li et al 2013, a trend that will likely continue under 21st century warming (Power et al 2013, Cai et al 2014, supporting continued research on the role of teleconnection patterns as drivers of terrestrial ecosystem processes.
Recent research efforts have identified important emerging links between vegetation activity and oceanatmosphere teleconnections at regional , Garcia and Townsend 2016), continental (de Beurs and Henebry 2008, Bastos et al 2016, Kim et al 2017, and global scales (Li et al 2012, Bastos et al 2013. While these efforts provide compelling evidence regarding the influence of teleconnections on vegetated ecosystems, critical knowledge gaps remain. For instance, many previous efforts considered only a single teleconnection seasonof-influence (usually prior winter) (e.g. Bastos et al 2016 and/or did not explicitly examine interactions among teleconnection patterns (e.g. . Given that teleconnection patterns are all embedded within a complex and interactive atmospheric system (Cayan et al 1998, improving knowledge of teleconnection interactions remains one of the key uncertainties in understanding how atmospheric circulation affects terrestrial vegetation . Our objective in this research is to improve understanding of atmospheric teleconnection influence on land surface phenology across North America, with explicit examinations of the importance of teleconnection interactions and the predictability of land surface phenology based on teleconnection indices. We analyzed the coupling of remotely sensed estimates of North American land surface phenology with the atmospheric components of five prominent teleconnection patterns. First, we used the Southern Oscillation Index (SOI), defined as a pressure difference between Tahiti and Darwin, Australia (supplementary figure S1 available at stacks.iop.org/ERL/13/034029/ mmedia), as an indicator of the atmospheric component of the coupled ocean-atmosphere ENSO system. Positive phases indicate a large pressure gradient across the equatorial Pacific (La Niña-type conditions) while negative phases indicate a weak pressure gradient (El Niño-type conditions). ENSO affects the position and intensity of Pacific storm tracks , with downstream impacts on both North American climate , Lee 2017) and vegetation activity (Bastos et al 2013. We also examined the two leading modes of north Pacific atmospheric variability: the Pacific-North American (PNA) and West Pacific (WP) patterns (figure S1). Positive PNA indicates ridging over western North America, resulting in higher temperatures and lower winter moisture delivery to western North America (Leathers et al 1991), while negative PNA indicates more zonal flow (Wallace andGutzler 1981, Barnston andLivezey 1987). Strong WP anomalies affect the position and intensity of the Pacific/East Asian jet stream (NOAA Climate Prediction Center 2012b), and changes in this system influence winter temperature and precipitation throughout much of North America (Linkin and Nigam 2008). Finally, we examined the two leading modes of atmospheric variability over the north Atlantic: the North Atlantic Oscillation (NAO) and East Atlantic (EA) pattern (figure S1). The NAO reflects a north-south pressure dipole between southern Greenland and the Azores (Barnston and Livezey 1987). Negative NAO phases typically lead to below normal temperatures in the eastern US (Osborn 2011). The EA, which is strongly associated with changes in the North Atlantic jet stream (Woollings et al 2010), reflects a northwest-southeast pressure dipole with centers displaced southeast of the NAO (Barnston and Livezey 1987).
Our goals in this research were to examine five key questions regarding the influence of these teleconnection patterns on North American land surface phenology: (1) Which of the teleconnections are most influential? (2) When do they exert the most influence?
(3) Where is phenology most responsive to atmospheric teleconnections? (4) Is there evidence of interactions between different teleconnection patterns? (5) How predictable is vegetation activity based on leading teleconnection patterns? Answering these questions could both improve short-term forecasting of ecosystem processes and illuminate potential long-term impacts of changes in the mean or variability of these circulation systems.

Teleconnection indices
We obtained monthly indices for all five teleconnections from the National Oceanic and Atmospheric Administration. Monthly SOI was calculated from standardized sea level pressure differences between Tahiti and Darwin, Australia (NOAA National Centers for Environmental Information 2017). The four northern hemisphere teleconnection indices represent leading modes of variability from a rotated principal component analysis (RPCA) of standardized monthly 500 mb geopotential height anomalies, with the RPCA procedure conducted separately for each month of the year to account for seasonal variation in the spatial patterns of atmospheric circulation (Barnston and Livezey 1987, NOAA Climate Prediction Center 2012a). We composited each index into 3 month running means, with subsequent analyses performed using all 3 month composites ranging from September-November (SON) prior to the northern hemisphere growing season through July-September (JAS) during the growing season.

Land surface phenology estimates
We developed four land surface phenology metrics covering 1982-2015 from the twice-monthly 3rd generation Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index (NDVI3g) (Pinzon and Tucker 2014): timing of the start and end of the growing season (SOS and EOS, respectively), total growing season length (LOS; days between SOS and EOS), and peak growing season NDVI (NDVI max ). We discarded all NDVI3g observations that were flagged as snow-or cloud-covered (Pinzon and Tucker 2014), as well as pixels that were composed of 50% or more of unvegetated land (barren or water) or non-natural vegetation (cropland and urban areas) based on a 1 km land cover classification (figure S2) (Hansen et al 2000). To further reduce noise and aid in visual interpretation, we spatially aggregated all NDVI3g time-series by taking the median value within each of the Level 3 EPA ecoregions of North America (Omernik 1987). We excluded any ecoregions with mean annual NDVI less than 0.1 or with no evidence of an annual phenological cycle in wavelet spectra (Torrence and Compo 1998). We also excluded any ecoregions with fewer than 25% of pixels remaining following the land cover filtering described above. For the remaining ecoregions, we used outlier and modified best index slope extraction filters to remove potentially unreliable NDVI3g observations and then fit a double-logistic curve to the annual filtered NDVI3g of each ecoregion (figure S3) (Hwang et al 2011a, 2011b. We identified phenological timing (SOS and EOS) using local half-amplitude thresholds and defined NDVI max as the annual maximum of the fitted doublelogistic curve (figure S3) (White et al 1997. Since half-amplitude thresholds would be affected both by changes in the initial onset of greenness as well as the timing and magnitude of NDVI max , we also performed additional analyses with a 10%-amplitude threshold, which would be mostly sensitive to the initial onset of greenness and minimally sensitive to the timing and magnitude of NDVI max . Of the 182 North American ecoregions, twelve were excluded due to lack of natural vegetation cover (figure S2), and an additional three were excluded due to no significant power at annual frequencies. Another eight ecoregions, mostly at very high latitudes, were excluded due to low mean annual NDVI. In total, the excluded ecoregions were mostly located in low productivity regions (the Arctic and the Southwest), predominantly agricultural regions of the Midwest, and regions with little seasonal variation of NDVI (e.g. southern Florida) ( figure S4). For the remaining ecoregions, the phenology of each individual year cannot always be resolved, particularly in dry regions where vegetation is sparse (often with short, rainfall pulse-driven, and/or bimodal growing seasons), and in densely-vegetated regions where NDVI is consistently high year-round (figure S4). We therefore restricted subsequent analyses to only include time-series of phenological metrics with at least 30 annual retrievals out of the 34 year study period. This results in a total of 129 ecoregions with SOS estimates (figure S4(a)), 115 with EOS estimates (figure S4(b)), and 149 with NDVI max estimates (figure S4(c)).

Statistical analyses
SOI, NAO, and EA had significant trends (p < 0.05) during the 1981-2015 period (figure S5), so relationships between teleconnections and land surface phenology could reflect common trends rather than explanatory power at interannual time-scales (Gonsamo et al 2016). We therefore detrended the monthly teleconnection indices during their period of overlap with the satellite data (1981-2015; figure S5) and formed 3 month seasonal composites from the monthly residuals. We also detrended the annual ecoregionscale phenology metrics using linear models, with the residuals used for all subsequent analyses.
We identified the dominant teleconnection seasonof-influence using seasonal correlation analyses, in which the detrended phenological metrics from each ecoregion were correlated against each of the 3 month teleconnection composites using Pearson's correlation coefficient (p < 0.05). Since multiple hypothesis testing on a large number of ecoregions can result in spurious rejections of the null hypothesis (i.e. false discovery), we also conducted 'field' significance tests which control the false discovery rate (Wilks 2006(Wilks , 2016. For a given teleconnection composite, we conducted n hypothesis tests, where n is the number of ecoregions tested, resulting in a spatial field of correlation values and their associated p-values. We compared the distribution of p-values from the ecoregion-level hypothesis tests to the distribution that would be expected if all null hypotheses were true, with the overall spatial correlation field considered significant if any local null hypothesis was rejected when controlling the false discovery rate at the 0.05 level (Wilks 2006(Wilks , 2016. We used a set of ordinary least squares multiple regression models to assess both the predictability of land surface phenology based on teleconnection patterns and the presence of interactions between pairs of teleconnections. For each ecoregion, we defined a set of additive (equation 1) and interaction (equation 2) regression models between pairs of 3 month teleconnection composites: Here,̂represents modeled phenology in year t, and x and z represent a pair of teleconnection indices used to predict̂. In equation 1, the effects of the teleconnections are additive: the phase of one does not modify the effect of the other. In equation 2, we include an interaction term ( 3 x t z ) which allows the effect of one teleconnection to be modulated by the phase of the other.
We assessed model skill and the importance of teleconnection interactions using both (1) partial Fstatistics to test for significance of the interaction term ( 3 ) and (2) leave-one-out cross-validation of the linear models, in which annual phenology estimates were successively withheld from model fitting. We evaluated model predictions against these withheld observations using the root mean squared error (RMSE) and the Nash-Sutcliffe model efficiency (E) (Nash and Sutcliffe 1970), defined as: where n is the number of annual observations, y is the observed phenology in year t, and̄is the mean phenology during the study period . E can vary from −inf to 1, with values greater than 0 indicating that the model is a more effective predictor of y than the observed mean value (̄) and a value of 1 indicating that the model is a perfect predictor of y ; models with E ≤ 0 have no predictive skill. In addition to the F-tests of 3 , we consider there to be evidence of an interaction between a given pair of teleconnections if E > 0 for the interaction model and if the RMSE of the interaction model is lower than the RMSE of the additive model.

Seasonality and spatial patterns of relationships between teleconnections and phenology
The three Pacific teleconnection patterns were each significantly correlated with SOS timing during at least one season (figure 1 and S6-S7), even when controlling for false discovery rate. Prior autumn-early spring conditions in the tropical Pacific and winter-spring conditions in the north Pacific exerted a particularly strong influence on SOS timing ( figure 1(a)). By comparison, the north Atlantic circulation patterns (NAO and EA) were not significantly related to SOS timing in most ecoregions (fewer than 15% for most temporal composites), and the false discovery tests suggest that these correlations were likely spurious (figures 1(d) and (g)). Few of the teleconnection indices were strongly correlated with EOS timing or with NDVI max (figure 1). We therefore focus solely on teleconnection impacts on SOS timing. SOI composites ranging from prior year SON through current year FMA were each significantly correlated (p < 0.05) with SOS timing for roughly 15%-21% of North American ecoregions (figures 1(a) and S6), and the overall correlation field was significant at the 95% level even when controlling for false discovery rate. Correlations between SOI and SOS were positive throughout the Pacific Northwest (figure S7), with effect sizes ranging from a couple days to more than a week of SOS delay per standard deviation increase in SOI (figure S8), likely due to higher temperatures in the Pacific Northwest following El Niño-type events (figure S9). The relationship between SOI and SOS timing resulted from combined effects on both the initial onset of greenness (figure S10(a)) and the timing of peak greenness (figure S10(c)). However, the earlier SOS following El Niño-type phases was partially canceled out by earlier EOS, as suggested by fewer significant correlations between SOI and LOS (figure S11(a)). The spatial patterns and sign of the correlations were similar for both prior autumn (OND) and late winter/early spring (FMA) composites of the SOI (figures S7 and S8).
Winter to early spring composites of the two north Pacific teleconnection patterns (PNA and WP) were significantly correlated with the SOS timing in roughly 15%-22% of ecoregions (figures 1(j) and (m)), with strongest correlations to PNA DJF and WP MAM . Correlations between PNA DJF and SOS were negative in and around the Great Lakes and southern Canada (i.e. earlier growing seasons following positive PNA DJF phases), while correlations between WP MAM and SOS Figure 1. Percentage of North American ecoregions with significant correlations (p < 0.05) to a given 3 month seasonal composite of each teleconnection index. X-axis labels represent the ending month of each composite, with asterisks ( * ) representing months in the prior year. For example, D * would represent a seasonal composite of October through December of the year prior to growth. Dark bars are shown where at least one ecoregion was significantly correlated with a given index when controlling the 'false discovery rate' at the 0.05 level (Wilks 2006(Wilks , 2016 were negative across much of the boreal region (figure S7). Effect sizes in these regions ranged from a few days to nearly one week earlier SOS per standard deviation increase in PNA DJF and WP MAM (figure S8), likely reflecting higher temperatures following positive phases of these north Pacific teleconnections ( figure  S9). The relationships between SOS timing and these north Pacific teleconnections resulted primarily from effects on the initial onset of greenness and less from the timing of peak greenness (figure S10). For WP MAM in particular, the earlier SOS following positive phases generally resulted in longer growing seasons (figure S11(d)).

Importance of teleconnection interactions
Based on results from figure 1, we further examined three Pacific teleconnection composites (prior SOI OND , PNA DJF , and WP MAM ) as well as their interactions with two Atlantic teleconnection composites (EA MAM and NAO JJA ). Importantly for subsequent regression analyses, these five teleconnection patterns were relatively independent of each other, with less than 10% shared variance between most pairs (table S1). The main exceptions, as suggested by previous research (e.g. Cayan et al 1998), were 11% shared variance between SOI OND and WP MAM and 28% between SOI OND and PNA DJF (table S1). Likewise, the detrended composites used subsequently in this study exhibited low autocorrelation from year-to-year (figure S12).
In addition to strong relationships between SOS timing and individual Pacific teleconnections, we found evidence of important interactions among teleconnection patterns (figure 2). The most prominent interaction effects were between winter-spring Pacific  patterns and summer NAO, suggesting that the effects of these teleconnections on ecosystem function may be greater or lesser than the sums of their constituent parts. Notably, while NAO was weakly correlated with SOS on its own, we found evidence of interaction effects between NAO JJA and Pacific teleconnections (figures 2 and S13). Based on the additive vs interaction model tests (figures 2(a) and S13), approximately 13% of the ecoregions examined in this study had evidence of an important negative interaction between NAO JJA and SOI OND (primarily in western and central Canada), and 17.1% had important interactions between NAO JJA and PNA DJF (primarily in eastern and central Canada). Partial F-tests yielded similar results ( figure 2(b)), particularly with a significant PNA-NAO interaction in 15.5% of North American ecoregions. These interactions were reduced when using a 10% NDVI3g threshold instead of a half-amplitude threshold to define SOS timing (figure S14), suggesting that the interactive effects of NAO JJA may be more closely related to NDVI max timing than to the initial onset of greenness, which would be expected since the JJA temporal composite occurs after the initial onset of greenness in these systems (figure S15).
SOS anomalies during extreme SOI OND , PNA DJF , and NAO JJA phase combinations (defined as the lower and upper terciles of each index) further revealed the nature of these interactions (figures 3-4). As suggested by the sign of 3 (figures 2 and 5(a) and supplementary text S1), the effects of SOI OND on SOS were magnified   ure 3(a)) while negative SOI OND was associated with earlier SOS (figure 3(c)). However, these patterns were amplified during negative NAO JJA phases (figures 3(b) and (d)). Dividing the SOI OND and NAO JJA into phases based on indices ±0 in order to increase the sample size within each phase showed similar patterns of stronger SOI OND effects during negative NAO JJA phases (figure S16).
This interaction effect was even stronger between PNA DJF and NAO JJA (figures 2, 4 and 5(b)). SOS anomalies were generally weaker with respect to PNA DJF phase when NAO JJA was positive (figures 4(a) and (c)), but very strong when NAO JJA was negative (figures 4(b) and (d)), possibly resulting from changes in the length of the green-up period related to NAO phase. Growing seasons started substantially earlier than average (by 1-2 weeks) across many boreal ecoregions following PNA+ winters and NAO− summers ( figure 4(b)), and much later than average (also by roughly 1-2 weeks) following PNA− winters and NAO− summers (figure 4(d)). As with the SOI, dividing NAO and PNA into phases based on indices ±0 resulted in lower overall anomalies, but with the same pattern of stronger PNA effects when NAO is negative (figure S17).
Skill of predictive models SOS timing could be partly predicted (0.2 < E < 0.45) in many boreal ecoregions with both additive and interaction models (figures 5(c)-(d), S18 and S19) based on pairs of atmospheric teleconnection patterns. Models that included autumn-to-winter SOI performed well in the Pacific Northwest, particularly in British Columbia ( figure 5(c)), while models that included WP MAM performed well across much of southern Canada (figures S18 and S19). In many cases, these models improved when including interaction terms in the models. Models that included PNA DJF , for example, performed better when the PNA:NAO interaction was included (figures 5(d), S18, and S19). For both the SOI:NAO and PNA:NAO models, there was a clear northeastsouthwest gradient in the RMSE (figures 5(e)-(f)), with RMSE of a week or less in many east-central Canadian and US ecoregions and up to a month in the drier, more variable ecoregions of the western US

Discussion
In this study, we show that start of growing season (SOS) timing can be predicted based on knowledge of atmospheric conditions over the Pacific and Atlantic Oceans, in some cases with several months of leadtime. Models that included prior autumn-winter SOI, for instance, were effective over much of the Pacific Northwest, while models that included prior winterspring WP or PNA were predictive of SOS timing throughout many boreal and Great Lakes ecoregions. Though atmospheric conditions over the north Atlantic were not important predictors of SOS on their own, we found relatively strong and widespread interaction effects of Atlantic and Pacific patterns on SOS, particularly between summer NAO and autumn to winter SOI and PNA. While models based solely on Pacific teleconnections performed reasonably well (0.2 < E < 0.45) across a large swath of the boreal region, including interaction effects with NAO JJA improved model performance, though including this interaction would also eliminate advanced forecasting since NAO JJA includes information from summer months.
The mechanism behind the NAO JJA interactions with Pacific patterns remains unknown. Unlike many teleconnections, the NAO is present in all seasons of the year (Barnston and Livezey 1987) and has been shown to impact summer hydroclimate in North America (Ruiz-Barradas andNigam 2005, Coleman andBudikova 2013). However, in much of this region, SOS timing (as defined with the half-amplitude threshold) typically occurs prior to the summer (JJA) NAO period (figures S15(a-c)). We thus suspect that the interaction effects observed here may reflect effects of the summer NAO on the timing and magnitude of NDVI max , which typically occurs around late August in the boreal region ( figure S15(d)). Any shifts in NDVI max timing would affect the length of the green-up period, which would therefore influence SOS timing defined based on a half-amplitude threshold. The relative lack of NAO JJA interaction effects in models of SOS timing based on a 10% NDVI3g threshold (which would be influenced primarily by the initial onset of greenness) suggests that this is a plausible explanation (figure S14).
The relationships between SOS timing and the winter-spring Pacific teleconnection patterns can likely be explained by coherent, large-scale influence of these patterns on early growing season temperatures (figure S9). Temperate and boreal spring phenology are strongly linked to temperature (Richardson et al 2013, Piao et al 2015, and remotely sensed NDVI observations are well-suited to capturing early season phenology in these regions of North America (White et al 2009). These factors likely explain why we find higher predictability of SOS timing than EOS timing. Autumn greenness and leaf area tend to be decoupled from canopy chemistry and physiological activity (Bauerle et al 2012, Croft et al 2014, so NDVI estimates of EOS tend to be several weeks later than the actual date of photosynthetic cessation . The relatively low teleconnectionbased predictability of EOS found in this study could therefore be due to a combination of three factors: (1) low quality retrieval of EOS timing from NDVI3g, (2) weak relationships between teleconnection patterns and late-season surface climate, and/or (3) complex and variable responses of EOS to temperature, drought, or photoperiod (Bauerle et al 2012, Richardson et al 2013, Xie et al 2015. Most of the teleconnection composites examined here were not significantly correlated with North American summer or autumn temperature, so (2) is certainly plausible, but the relative contributions of the above explanations to weak EOS-teleconnection relationships remain poorly defined. It may be possible to improve EOS retrievals using more physiologicallybased remote sensing products, such as vegetation solar-induced fluorescence (Joiner et al 2014, which could further clarify the drivers of EOS timing. Accurate phenological forecasts at seasonal to interannual time-scales are important for land management, agriculture, and even for human health (e.g. through pollen and disease vector abundance) (Morisette et al 2009). Our results suggest that assimilating teleconnection indices into current short-term phenology forecasting systems based on autoregressive properties of vegetation index time-series (e.g. White and Nemani 2006) could represent a possible mechanism to improve seasonal forecasts of vegetation activity. Pacific teleconnection patterns in particular could improve forecasts of SOS timing across much of the boreal zone. Accounting for interaction effects with the Atlantic patterns could further improve SOS predictions, but with less lead-time since the strongest interactions were observed with later Atlantic teleconnection composites. Further research could examine influences of lower frequency oceanic oscillations like the Pacific Decadal Oscillation or the Atlantic Multidecadal Oscillation on land surface phenology, including whether these oceanic oscillations interact with higher frequency atmospheric circulation patterns. Oceanic teleconnection patterns tend to have longer persistence than atmospheric circulation (McCabe et al 2012), suggesting that they may provide even longer forecast lead-times than the teleconnection patterns examined here.
Examining these lower frequency, (multi-)decadal oscillations would likely require longer-term estimates of vegetation activity than those provided by remote sensing, which could include tree-ring records or land surface model ensembles like TRENDY (Sitch et al 2015) and MsTMIP (Huntzinger et al 2013). These datasets could also be used to address two limitations of our research. First, the remotely sensed data used here often fail to capture phenology and productivity dynamics across dry regions like the North American Southwest (White et al 2009). Growing seasons in dry regions tend to be driven by pulses of rainfall, leading to temporally-variable and often bimodal growing seasons which, in combination with the high soil background reflectance, complicates retrieval of vegetation information from NDVI (Biederman et al 2017). The lack of teleconnection influence on Southwestern vegetation found in this study may therefore reflect inherent limitations of the remote sensing data in this region. Second, while we focus solely on vegetation dynamics related to phenology and the annual uptake of CO 2 , model ensembles and CO 2 inversions could be used to examine the impacts of atmospheric circulation on total carbon flux, including both photosynthetic uptake and respiratory loss of CO 2 . Quantification of teleconnection influence on both the uptake and respiration sides of the carbon balance could add key insight into the mechanisms driving fluctuations in the atmospheric CO 2 growth rate, a current area of high uncertainty (Keenan et al 2016, Ballantyne et al 2017.

Conclusions
Our main objectives in this research were to determine the role of teleconnection interactions in regulating land surface phenology and to explicitly quantify the predictability of land surface phenology based solely on knowledge of teleconnection patterns. Regarding these objectives, there are five key findings of this research: 1. The tropical Pacific (SOI) and north Pacific (PNA and WP) patterns were more closely related to North American spring phenology than the north Atlantic indices (NAO and EA).
2. Pacific teleconnection patterns provided significant information on start of growing season timing several months in advance of the growing season, with prior autumn-winter conditions over the tropical Pacific and winter-spring conditions over the north Pacific being strongly correlated with start of growing season timing in many ecoregions.
3. Relationships between teleconnection patterns and vegetation activity were generally strongest across the energy-limited boreal and Great Lakes regions.
4. We found widespread evidence of interaction effects among teleconnection patterns, particularly between circulation patterns over the Pacific and those over the north Atlantic, even though the Atlantic patterns did not explain much variance in start of growing season timing on their own. 5. Multiple regression models based on pairs of teleconnection indices, particularly those that included at least one index of winter-spring atmospheric conditions over the Pacific, were effective predictors for start of growing season timing across much of the boreal zone, suggesting potential for these indicators to forecast growing season dynamics several months in advance.
Overall, our findings suggest that knowledge of atmospheric circulation patterns significantly increases the predictability of North American spring phenology. A more robust consideration of the complexity of the atmospheric circulation system, including interactions across multiple ocean basins, is an important step towards accurate forecasts of vegetation activity.