Intensification and poleward expansion of upwelling-favourable winds have been predicted as a response to anthropogenic global climate change and have recently been documented in most of the world. To identify how these processes are impacting nearshore oceanographic habitats and, especially, long-term trends of primary productivity in the Humboldt Upwelling Ecosystem (HUE), we * Documenting long-term trends in (EBUSs) is of great importance not only because a large fraction of the human population lives near these shores, but also because these are the most productive marine regions in terms of phytoplankton biomass and fish stocks (Cushing 1971;Cury et al., 1998;Chavez and Messie, 2009). Indeed, equator-ward winds along these shores induce the upwelling of subsurface cold water masses rich in nutrients, which fuel phytoplankton growth and, through the trophic web, provide the necessary energy to sustain some of the largest industrial and artisanal fisheries of the world (Mann and Lazier, 2006;Chenillat et al., 2013;Chavez and Messie, 2009;Salas et al., 2011). Thus, variability of climatic drivers on surface water productivity within EBUS regions impinges directly on our ability to adapt to climate change and on the sustainability of major industrial and local fisheries, with the large direct and indirect socio-economic impacts they generate (Pauli and Christensen, 1995). Consequently, several recent studies have examined long-term trends in wind forcing and physical oceanographic conditions in most EBUSs and, naturally, our understanding of long-term climate effects in these regions of the worlds' oceans has improved considerably. The effects that these long-term hydrographic changes are having on productivity of the coastal oceans are now a matter of intensive research. However, to understand these temporal trends in primary productivity, studies simultaneously integrating large scale climatic forcing, wind dynamics and sea surface temperature variability are still lacking for the most important EBUSs. On a global scale, the poleward expansion of the atmospheric Hadley circulation cells, due to the weakening of the thermal differences between the poles and the equator, is known to have already caused the poleward displacement of both the subtropical anticyclones and the sub-polar westerlies belts (Previdi and Liepert, 2007;Nguyen et al., 2013). The result has been a poleward expansion and intensification of upwelling favourable winds along most EBUS over the past two decades (McGregor et al., 2007;Lima and Wethey, 2012;Sydeman et al., 2014). The trend is particularly evident in the Southern Pacific, along the Humboldt Upwelling Ecosystem (HUE), where the poleward movement of the Southwestern Pacific Anticyclone (SPA) has displaced the westerlies belt closer to Antarctica (Fan et al, 2014;Ancapichun and Garces-Vargas 2015). Recent studies have shown that the seasonal latitudinal migration of the SPA has moved poleward to sit around 36 deg S during the critical austral spring-summer months that fuel phytoplankton productivity in the coastal section of the south eastern Pacific ocean (Schneider et al., 2017;Aguirre et al., 2018). The displacement has been captured by one of the main modes of climatic variability in the region, the Pacific Annular Mode (SAM), which has shown a clear positive trend in the last decades, thus pointing to low atmospheric pressures over Antarctica compared to those at subtropical latitudes (Marshall, 2003;Wang and Cai, 2013). It has also driven a shift in atmospheric conditions along a broad region of the central coast of Chile, between ca. 30-35 deg S, where spring winds have persistently decreased over the past decades (Aguirre et al., 2018), and a poleward region where winds have intensified, rainfall has decreased and surface water cooling has been observed throughout the mixed layer, especially since 2007 (Schneider et al., 2017;Jacob et al., 2018;Narvaez et al., 2019). The latitudinal position of the SPA and its regular seasonal migration clearly modulate long-term trends in upwelling phenology and atmospheric climatic conditions (Ancapichun and Garces-Vargas, 2015;Weller 2015, Jacob et al., 2018), and impose geographic discontinuities in upwelling regimes, which in the HUE occur around 30 deg S (Strub et al., 1998;Hormazabal et al, 2001;Navarrete et al., 2005;Tapia et al., 2014). Such discontinuities in upwelling regimes are known to have farreaching consequences on pelagic and benthic ecosystems, including population dynamics, larval recruitment, adult abundances, the role of species interactions, and genetic and functional structure of benthic communities (Navarrete et al., 2005;Wieters et al., 2009;Tapia et al., 2014;Haye et al., 2014). How climate forcing and the progressive poleward displacement of the SPA will in the long-term affect the biogeographic discontinuity at 30 deg S remains unclear. Besides the intensification and poleward migration of subtropical anticyclones, in 1990 Bakun (Bakun 1990) predicted that the greenhouse effect would warm the oceans slower than adjacent land masses, thus the land-sea pressure differential will become reinforced, strengthening landward breezes which, due to Coriolis, would increase equatorward winds and coastal upwelling at all EBUSs. The phenomenon is known as the Bakun effect (Bakun 1990;Bakun et al., 2010), and empirical evidence and model results supporting such response to anthropogenic warming have been controversial (Demarcq 2009;Pardo et al., 2011). Recent meta-analyses of observations and models show that significant intensification of winds, attributable to Bakun effect, has occurred in all EBUSs but the Canary Current system, where winds have weakened (Sydeman et al., 2014). This result is consistent with several in-situ observations (Garcia-Reyes and Largier, 2010;Weller 2015) and model predictions for the XXI century (Wang et al. 2015). But other studies, using different databases, have found no evidence for a general intensification of upwelling along EBUSs (Varela et al., 2015). Contradictory results obtained by different studies are due, at least in part, to the fact that long-term responses within EBUSs are not homogeneous, but trends may change across latitudes (Varela et al., 2015). Thus, the regional averages examined to date may provide an incomplete picture of long-term trends in the highly dynamic and spatially variable coastal ocean that characterize all EBUS. To what extent the positive or negative trends in equatorward winds, the driving force of upwelling, have led to changes in primary productivity of coastal waters is now being intensively investigated (Bakun et al., 2015;Gomez-Letona et al., 2017).
The effect of altered upwelling winds on surface primary productivity can have contrasting effects over coastal and offshore domains and analyses must therefore take this into account. For instance, stronger coastal upwelling generates turbulence and mixing alongshore which, together with enhanced seaward transport by the Ekman drift, may lead to higher offshore phytoplankton concentrations at the expense of a decrease in nearshore primary productivity (Lackhar and Gruber, 2012;Anabalon et al., 2016). In contrast, a scenario with higher phytoplankton biomass nearshore induced by enhanced upwelling-driven nutrient inputs is also possible (Bakun et al., 2010). Moreover, the influence of upwelling intensification or weakening on phytoplankton biomass may depend on the ecophysiological requirements of different phytoplankton groups (Smayda, 2000) or changes in other nutrient sources, such as riverine inputs (Masotti et al., 2018;Jacob et al., 2018). In this study, we aim to quantify and disentangle the role of SPA dynamics and Bakun effect on upwelling spatiotemporal fluctuations at the Humboldt Upwelling System. To this end, time series of sea-land thermal differentials, sea level atmospheric pressure, meridional wind stress, satellite sea surface temperature (SST) and chlorophyll-a concentration (Chl-a) with the highest possible spatiotemporal resolution for the region, together with in situ measurements of nearshore Chl-a, were inspected along central Chile between 26o-36 deg S in search of interannual and seasonal trends. With this approach we were able to evaluate the strength of such trends, their spatial distribution across the region and their correspondence with the displacement of the SPA and Bakun effect. The direction and spatial structure of the temporal trends measured in all these variables will allow us further comparisons with other models and observations developed both within the Humboldt Upwelling System and in other EBUSs. Processed, science-quality satellite data for SST and Chl-a were retrieved from the (MODIS) onboard the Aqua spacecraft 1 . In particular, 8day composites with 4 km of pixel resolution were retrieved at different distances from the shore (20,40,60,80,100,200,300,400 and 500 km) between 26o-36 deg S off the Chilean coast from May 2003 to December 2015. To reduce the spatial variability in the datasets and the gaps along the time series, latitudinal moving averages from North to South every 4 km using pixel aggregations (superpixels) with a length of 12 km alongshore  20 km cross shore (3  5 pixels) were calculated and used in further analyses. The ERA-Interim weather reanalysis model from the European Center for Medium-Range Weather Forecast ( was used to obtain meridional wind velocities at a height of 10 m (V ). Input data for wind data prediction and reanalysis is provided by Meteosat satellites (Dee et al., 2011). Model estimates were downloaded as monthly means for pixels of 0.75 deg from 26.25o S to 36 deg S and from 73o W to 69.75 deg W. Then, meridional wind stress, the component which drives upwelling along the Chilean coast, was calculated asas it has been previously used to characterize Ekman transport dynamics in the region (Aguirre et al. 2014). It was computed as follows: where W stress stands for wind stress,  air is air density (1.225 kg m -3 ), CD is the drag coefficient (0.0014) and V is meridional wind speed. Positive and negative values represent northwards and southwards wind forcing, respectively. In order to infer SPA position over time, from the ERA Interim model monthly sea level atmospheric pressures at pixels with a size of 0.75 deg were retrieved from 2003 to 2015 in between 76o-96 deg W and 25o-40 deg S. Then, a third-order polynomial fit was applied to sea level pressures as a function of latitude following the method of Minetti et al. (2009). The SPA position was calculated for the entire time series as the latitude where the maximum of the polynomial curve (first derivative zero) is reached. To evaluate Bakun effect across the study region over the timespan encompassed by our study, monthly ERA Interim modelled surface atmosphere ('skin') temperatures at pixels with a size of 0.75 deg were retrieved from 2003 to 2015 from the first inland pixel in between 26o-36 deg S. Then, the monthly sea-land thermal differences (in deg C) was obtained by subtracting monthly average SST values at the 20 km superpixel and over 0.75 deg latitudinal bands, from the monthly land skin temperatures. Chlorophyll-a in surface waters: in situ and satellite data Surface Chl-a data were retrieved from the (MODIS) at the Aqua spacecraft using the same procedure described above for SST. To validate satellite data and examine in-situ time trends, we used the time series of daily, Chlorophyll-a extracted at the Estacion Costera de Investigaciones Marinas (ECIM) located at 33.49 deg S . Three 250 ml water samples were collected daily at depths between 20-40 cm from 2003 to 2012 and filtered onto 0.2 um filters. Then these were prefiltered with a 300m mesh and 100ml were filtered through a GF/F Whatman glass fibre filter (Wieters et al., 2003). The Chl -a was extracted by adding 5 ml of 90% acetone during 24 h at -20 deg C to the GF/f filters (Yentsch and Menzel, 1963). Concentrations were then measured with a 10-AU Turner Designs fluorometer and used to obtain daily averages for a time series. (GAMMs, Chen, 2000;Wood, 2006) were used to decompose time series of wind stress estimates, SPA positions, sea-land thermal differences, SST, Chl-a, and in situ Chl-a at ECIM, into a seasonal and a long-term trend component. Like standard Generalized Additive Models, GAMM fits a penalized cubic spline smoothing functions through backfitting procedures, but it allows to account for temporal autocorrelation and lag effects in the error structure through the inclusion of an autoregressive term, which is considered a random effect in the GAMM. The model is robust to data gaps and irregular time intervals between measurements along the time series (Simpson, 2018).
These analyses were implemented with the "mgcv" package (Wood, 2006) for R statistical computing version 3.4.0 (R Core Team 2019) and have been used before to analyse time series of hydrographical and biological variables (Harding et al., 2016). A linear fit was applied to the long-term trend component to obtain a time integrated slope or rate of change for wind stress, SST and Chl-a. In order to infer spatial patterns in temporal trends these slopes were calculated for each smoothed superpixel of 12  20 km in the case of MODIS data and for each 0.75 deg pixel in the case of ERA Interim wind stress estimates and represented in maps for the whole region. To examine the coherence of the long-term trends among seasons, especially during the period of maximum primary productivity in spring, linear trends for each season were performed separately for each wind stress pixel and MODIS superpixel. Seasonal means for each variable were obtained using the following periods: spring (September to November), summer (December to February), autumn (March to May) and winter (June to August). Seasonal slopes of these linear trends were also represented in maps for the whole region. Since previous studies have reported that poleward displacement of SPA has accelerated since 2007 (Schneider et al., 2017;Aguirre et al., 2018), mean seasonal sea level pressure maps for the periods 2003-2007 and 2008-2015 were compared in search of contrasting spatial patterns in the location of the SPA. We also calculated the difference in meridional wind stress, SST and Chl-a between the two periods for those coastal areas showing the most contrasting long term trends in Chl-a within the study domain. Chl-a and SST data were averaged for the first 100 km offshore at these coastal areas and the respective 0.75 deg pixels were selected to calculate the differentials of wind stress between the periods before and after 2007. These differentials were compared with those observed in previous studies further south at 36.5 deg S (Schneider et al., 2017;Aguirre et al., 2018). Time series of monthly SPA latitudinal position revealed strong seasonality, with the core of the anticyclone located at 29 deg S in austral winter (July and August) We observed highly significant positive long-term trends in sea-land thermal differences (i.e. Bakun effect), yet the temporal pattern showed a clear latitudinal discontinuity. North of about 30 deg S, a significant increase took place in both annual and spring sea-land pressure gradient, with a slope of about 0.08 deg C year -1 over the timespan considered. South of 30 deg S, annual pressure differential showed no trend and the slope remained close to 0, while spring pressure differential showed a non-significant positive trend at some latitudinal bands . On average, across the entire study region, long-term trends of meridional wind stress were mostly positive for the entire year at around 0.001 N m -2 year -1 , with slightly steeper trends north from 28 deg S and nearly null or even slightly negative trends (-0.0005 N m -2 year -1 ) close to shore between 30o-34 deg S. During the spring upwelling months, two contrasting regions in terms of trends in wind stress were clearly defined. Firstly, a region or domain north of 30 deg S, where marked increases in meridional wind stress took place at a rate of more than 0.002 N m -2 year -1 at 27 deg S. Secondly, a region or domain south of 30 deg S to 35 deg S, where meridional winds decreased at rates around 0.0015 N m -2 year -1 . Intensification of upwelling winds in the northern domain also took place during summer, while in autumn the increasing trend was only evident offshore (<100 km off the coast) in between 31o S to 34 deg S (S1). The response of SST yearly averages over the study period has been a coastal warming of around 0.05 deg C year -1 in the region south of ca. 29 deg S and down to 35 deg S , reaching maximum values of 0.1 deg C year -1 just north of Punta Lengua de Vaca (e.g. 30 deg S, Fig. 3B). Further north, around 27 deg S, yearly average temperatures have decreased at rates of 0.075 deg C year -1 . In austral spring months, long-term trends indicate a warming rate of 0.05 deg C year -1 along most of the domain, except for a coastal band north of about 28 deg S, where surface cooling has taken place . Similar cooling patterns were observed for austral winter months further offshore (>80 km offshore, S1). Long-term yearly trends in surface Chl-a emphasize the contrast between a northern domain, between 26 deg S and 30 deg S and a southern domain that extends to 35 deg S. Equatorward, yearly primary productivity has increased at rates of 0.4 mg m -3 year -1 , while in the poleward domain it has decreased at faster rates (-0.6 mg m -3 year -1 , Fig. 2C-D). A transition region of little or no temporal trends in overall Chl-a persisted around 30o-31 deg S . Spring months showed a broadly similar pattern marked by mesoscale spatial structure. Sharp long-term decreases in Chl-a with rates steeper than -1.2 mg m -3 year -1 were widespread in the southern domain, especially between 34o-35 deg S , while the northern domain showed increases of about 0.4 mg m -3 year -1 , for example between 29o-28 deg S, but other coastal areas show little or no change . Similar contrasting trends between the two latitudinal domains were observed for the rest of the seasons (S1). Daily time series of in situ Chl-a at ECIM were highly consistent with the pattern in primary productivity found at the southern domain from satellite data. Starting in 2008, the characteristic spring and fall blooms that were apparent in the time series between1999-2008, had all but disappeared , generating a marked long-term decrease in situ Chl-a concentration at this coastal site.
Although a negative linear function significantly fits the long term GAMM trend, with a slope of -0.22 mg m -3 year -1 (p-value <0.0001), the temporal trend at the site was more complex, with a clear shift to overall lower values after a sharp peak in austral spring 2007-2008 (Figs. 3I,K) (Figs. 3I,K). The seasonal component of the GAMM was also significant (p-value <0.0001) and showed a marked peak in October (spring) followed by a progressive decrease through summer and autumn months toward minimum winter values in July . Differentials in wind stress forcing, SST, and the consequent changes in Chla, between the period before and after 2007, which was the year that SPA poleward migration apparently accelerated (Schneider et al., 2017), provide only weak support for the coastal upwelling driven productivity model. The northern coastal region presented a marked increase of 0.0038 N m -2 in wind stress between these two periods, but a nil or even a slight rise of SST of about 0.005 deg C, yet a significant increase in surface Chl-a of 0.6 mg m -3 did occur . At the southern coastal domain, southerly wind stress dropped 0.0002 N m -2 between these periods, SST rose 0.25 deg C and Chl-a decreased 1.5 mg m -3 , thus closely following the responses expected from the upwelling-driven productivity model . Further south, at 36.5 deg S where a long-term monitoring site is located (Schneider et al., 2017), we did find an increase in meridional wind stress, as shown in previous studies, although the magnitude was smaller than 1e-4 N m -2 . Surface water temperature dropped 0.25 deg C between these two periods, but instead of increasing, Chl-a fell 0.75mg m -3 . Here we documented the existence of contrasting long-term trends in wind forcing and SST along different portions of the Humboldt Upwelling Ecosystem, and show how these climate-driven changes have translated into dissimilar domains of surface water productivity in the costal ocean. At the equatorward domain of the region, north of 30 deg S, phytoplankton productivity (Chl-a biomass) has increased in nearshore habitats since the early 2000's. At the poleward section, between 31o-35 deg S and over the same time period, phytoplankton productivity has significantly decreased over an extensive cross-shore region encompassing the entire continental shelf. The upwelling-driven surface productivity model provides a partial explanation for the observed trends in Chl-a productivity in the two domains of the region, evidencing the complexity of ecosystem responses to altered upwelling winds (Bakun et al., 2015). The long-term changes in wind forcing are in remarkable agreement with the temporal and latitudinal changes in two processes linked to anthropogenic global change in most EBUS, including HUE. First, we document an increased sea-land temperature gradient (i.e. Bakun effect) north of 30 deg S since the early 2000's that has led to increased upwelling winds in that portion of the region. Second, the well documented poleward displacement of the SPA beyond 35 deg S has increased spring winds south of that latitude but created a gap in between 31o-35 deg S where upwelling winds are now significantly weaker. The regular seasonal changes in the position of the SPA are known to drive upwelling variability along most of the HUE (Strub et al., 1998;Ancapichun and Garces-Vargas 2015). During the second half of the XX century, the core of the SPA in spring months was located at 33 deg S migrating farther south in summer and back up north to its autumn/winter position around 26 deg S . Yet over the past two decades, and in agreement with patterns observed for the high pressure systems of the North and South Atlantic (Kim et al., 2015;He et al., 2017;Reboita et al., 2019), the spring and, to a lesser extent, summer positions of SPA have been displaced poleward, reaching around 35 deg S and 37 deg S in spring and summer, respectively (Figs. 1A). Previous studies using combinations of ERA Interim reanalysis and GCM multimodel simulations (Aguirre et al., 2018), together with primary data from QuikScat and ASCAT scatterometers (Schneider et al., 2017), have also documented the increased SPA poleward spring displacement, lending support to our results. The displacement has led to the intensification of spring winds at the southern-end of HUE (Jacob et al., 2018;Narvaez et al., 2019), and a region of weaker spring winds in central Chile, between 30o-35 deg S. Decreased upwelling winds in this extensive region of the Chilean coast has caused a weak but significant trend to increasing SST, which is particularly apparent in front of the upwelling hot-spots of Matanzas (34.1 deg S), Punta Toro (33.8 deg S) and at the well-documented upwelling centre of Punta Lengua de Vaca (30.2 deg S). As expected from the upwelling-driven coastal productivity model (e.g. Botsford et al., 2003), weakening upwelling winds have led to a regional domain with decreased primary productivity for coastal waters over the past two decades, which is evidenced by long-term patterns of spring and mean annual Chl-a from satellite and in situ records. The region of decreased productivity extends well beyond the coastal zone, especially in spring months when surface waters with significantly reduced concentration of Chl-a are found as far as 100 km offshore . Such an extensive region of lower annual and spring Chl-a is bound to have effects on higher trophic levels, as the reduced primary productivity cascades up the pelagic and benthic food webs (Sydeman et al., 2006;Powell and Xu 2011). Regional displacements of coastal fisheries, top predators and sea birds to areas of higher productivity should already be taking place or should be observed in the near future. However, the decreasing Chl-a rates are not spatially homogeneous within the study region, exhibiting important and ecologically-significant meso-scale variation (Navarrete et al., 2005;Valdivia et al., 2015), with small nearshore sectors that have experienced almost no change (rates of less than -0.1 mg Chl-a / yr) to regions of very rapid decrease (-1 mg Chla / yr).
The long-term decreasing trends in Chl-a detected from satellite imagery few kilometres offshore (first pixel 20 km from shore) have also been observed in time series of in situ chlorophyll-a recorded in water samples collected from the intertidal zone at ECIM . Thus, the domain of decreased productivity includes nearshore waters and it should therefore affect the rich and heavily exploited intertidal and subtidal benthic ecosystems (Kefi et al., 2015, Perez-Matus et al., 2017). Equatorward of 30 deg S, we show that upwelling-favourable winds have actually shown significant increases, especially in spring months. Such wind intensification is consistent with the increased sea-land thermal differences observed from north of 30 deg S , thus we interpret the spatial pattern as a geographically constrained Bakun effect (Bakun, 1990). Note that upwelling strengthening in the northern domain cannot be explained by SPA dynamics, as in spring the SPA is located further south. This spatial decoupling between migration of pressure systems and Bakun effect has been observed in Monterey Bay, in the California Current System, where sea-land thermal diurnal breezes cause pronounced upwelling events regardless of the prevalent forcing exerted by the North Pacific High off California (Woodson et al., 2007). Increased winds in our northern region have not led to generally decreasing SST trends . Long term cooling in the northern region has apparently occurred only locally (between 27o-28 deg S), while many sectors have experienced either no change or even weak warming, probably due to the complex interaction between the upwelling dynamics, increased solar radiation (Echevin et al., 2012;Aguirre et al., 2018) and changes in other processes that modulate thermal conditions and which also change around this latitude (Tapia et al., 2014). For instance, stronger stratification by atmospheric warming may lead to upward transport of warm and already nutrient depleted waters to the surface unable to trigger conspicuous phytoplankton blooms. This scenario has been forecasted for the 21 st century within the HUE and off the Northwestern Iberian Peninsula, as increased stratification will lead to a shallower thermocline (Oyarzun and Brierley, 2018;Sousa et al., 2020). Even if meridional wind stress increases, a shallow thermocline may trap wind induced turbulence and currents, thus upwelled waters may ascend from shallower depths (Pringle 2007). It is worth noting that our results contrast with earlier analyses of long-term satellite trends in coastal SST over our study region, which showed a largely homogenous decreasing trend over the 1982-2010 period (Lima and Wethey, 2012). The difference can be reconciled with the shift in the SPA position observed from 2010 onwards. Overall, increases or decreases in upwelling winds and the position of the highpressure systems do not only alter offshore advection, but also many other processes that can affect productivity of the coastal ocean (Bakun et al., 2015). For instance, increased upwelling winds at 36.5 deg S (Station 18) due to SPA displacement, have produced sharp cooling of the water column, but Chl-a productivity has decreased . It has been hypothesized that lower precipitations due to droughts linked with the SPA intensification at that latitude, may be negatively affecting riverine outflow and reducing primary productivity (Jacob et al., 2018). Thus, the climatic drivers of upwelling intensity do not only alter offshore advective transport, but several other factors that also modulate ocean productivity, and riverine input is one of important one in several regions of EBUS (Thomas and Strubb 2001;Silva et al., 2009;Masotti et al., 2018). The driving mechanisms underlying the two contrasting domains in primary productivity of the coastal ocean along the HUE, as we propose here, can be associated to anthropogenic global warming. First, although the dynamics of the SPA are partially driven by the natural variability of the Pacific Decadal Oscillation (Ancapichun and Garces Vargas, 2015), anthropogenic warming appears to be the main cause of the poleward displacement of high EBUS pressure systems (Lee 1999;Garreaud and Falvey, 2009;Wang and Cai 2013;Fan et al., 2014;Aguirre et al., 2018). Second, Bakun effect is a response to the increased sea-land pressure differential produced by global warming, which should be intensified at lower latitudes (Bakun, 1990). We show that the Bakun effect does not change gradually with latitude but exhibits a sharp discontinuity at about 30 deg S. The sharp increase of the Bakun effect equatorward of 30 deg S is likely related to stronger land warming over the coastal section of the arid Andean slopes (Vargas et al., 2007). Such land warming may be weaker to the south of 30 deg S, where coastal lowlands lay in between the Pacific shores and the Andes. Further south and using data from station 18 at 36.5 deg S, complemented with wind data from satellite scatterometers, Schneider et al. (2017) showed an abrupt acceleration of SPA poleward displacement since 2007. Our analyses showed important differences in wind stress, SST and especially Chl-a between the beginning of the 2000's and the last decade of the time series at the 2 domains studied here and the southern end (station 18) analysed in previous studies . But SPA latitudinal position changed progressively over time and most of the time trends detected by the GAMMs all over the region were linear with no evidence of discontinuity. However, a clear exception was the longterm trend of in situ Chl-a at ECIM, which showed a marked decrease only after a conspicuous Chl-a peak in spring 2007 . This temporal pattern at ECIM was mirrored by the outflow of the Biobio river further south that also dropped after a large peak in 2007 (Jacob et al., 2018), suggesting large climaticchanges took place on that year. Cold SST anomalies nearshore during 2007 up to 29 deg S were thought to be driven by a change from El Nino to La Nina conditions or by a northward intrusion of sub-Antarctic waters (Aravena et al., 2014;Schneider et al., 2017). How these processes may interact with SPA dynamics to produce more or less pronounced shifts in oceanographic conditions all along the region remains unclear.
In summary, two main effects of anthropogenic global warming on upwelling systems seem to be responsible for two contrasting domains of long-term dynamics in primary productivity of the coast, with limits between these domains around 30o S to 31 deg S. To the north, Bakun's increased sea-land thermal contrast would enhance upwelling-favourable winds and coastal primary productivity. To the south of the latitudinal break, there is an ample region of decreasing upwelling winds that extends down to 34o-35 deg S due to the poleward displacement of SPA, where Chl-a concentrations have been drastically reduced. Thus, anthropogenic climatic trends may generate much more marked and conspicuous spatial heterogeneity in upwelling intensity and coastal water productivity than previously thought, which will affect entire coastal ecosystems by propagation through the food web. Our results highlight the urgent need for monitoring the coastal ocean to assess existing trends and forecast future changes in productivity within EBUS, and in this manner increase our capacity to adapt human activities to the cascading and far reaching effects of global warming. We thank research assistants and students, especially Ivan Albornoz, for helping in the collection of water samples at ECIM and Mirtala Parrague for managing the in-situ Chlorophyll-a dataset. The National Fund for Scientific and Technological Development, FONDECYT (Chile), supported with post-doctoral Conversations with Catalina Aguirre and colleagues at ECIM are much appreciated. Additionally, AO was SAN was We thank research assistants and students, especially Ivan Albornoz, for helping in the collection of water samples at ECIM and Mirtala Parrague for managing the in-situ Chlorophyll-a dataset. The National Fund for Scientific and Technological Development, FONDECYT (Chile), supported with post-doctoral Conversations with Catalina Aguirre and colleagues at ECIM are much appreciated. Additionally, AO was SAN was Preprint Preprint