Historical dynamics in ecosystem service bundles Delphine Renarda,1, Jeanine M. Rhemtullab, and Elena M. Bennettc a Department of Geography and Natural Resource Sciences, McGill University, Montreal, QC, Canada H3A 2T5; bDepartment of Forest and Conservation Sciences, University of British Columbia, Vancouver, BC, Canada V6T 1Z4; and cDepartment of Natural Resource Sciences and McGill School of Environment, McGill University, Ste. Anne de Bellevue, QC, Canada H9X 3V9

Edited by Monica G. Turner, University of Wisconsin-Madison, Madison, WI, and approved September 16, 2015 (received for review February 6, 2015)

|

|

|

ecosystem services historical ecology bundles ecosystem service interactions spatiotemporal analysis

|

M

anaging multiple ecosystem services (ES) simultaneously, including addressing trade-offs between services and preventing ecological surprises, is among the most pressing concerns of sustainability research (1–3). However, most ES research to date cannot truly address these critical challenges because it has focused primarily on quantifying and mapping the delivery of only a few services at a single point in time (4). In this study, we analyze nine ES at five-year intervals from 1971 to 2006 to show how a spatiotemporal approach can enhance our understanding of ES dynamics. The adoption of a historical perspective has made important contributions in other areas of ecology (5–8). For example, time has been revealed to be as important as space for understanding patterns of species richness and distribution (6, 9). Historical ecology has shed light on the persistent effects of human activity on landscapes (10) and ecosystem function (11–13). This field has also provided the temporal perspective needed to understand the underlying causes and rates of change in ecosystems as context for the future, including the likelihood of unexpected regime shifts (14–16), and the potential for conservation, restoration, and management of ecosystems (17–19). However, historical analyses have been largely absent from ES research thus far. The few studies presenting a historical approach mostly compare two snapshots in time (20) or quantify only a limited number of ES (21, 22), and rarely investigate interactions among multiple services through time (but see (23)). These shortcomings are of particular concern given that the demand for most services is increasing, making interactions such as synergies

www.pnas.org/cgi/doi/10.1073/pnas.1502565112

and trade-offs among ES more important to take into account. A better understanding of how multiple ES interact, how trade-offs and synergies emerge, and how interactions may shift through time as conditions change or respond to new drivers can help meet the challenge of managing multiple ES. In this study, we examined how multiple ES and the relationships among them have changed over time and across space. We extended the spatially explicit ES-bundle approach of RaudseppHearne et al. (24) with data spanning 35 y. We defined an ES bundle as a mix of positively correlated ES provided together in the same place and at the same time, even though they may not have causative relationships. Using primary data compiled from diverse sources, we quantified the delivery of nine ES (Table 1; four provisioning, two regulating, and three cultural services) at 5-y intervals from 1971 to 2006 and across 131 municipalities in the Montérégie, a suburbanizing and heavily agricultural region in southern Quebec in Canada (11,853 km2). Our multitemporal and spatially explicit approach provided an opportunity to (i) assess the extent of temporal and spatial variation in nine ES individually; (ii) identify ES bundles, examine their spatiotemporal dynamics and how they relate to environmental and socioeconomic characteristics of the study region; and (iii) determine changes in the relationships (trade-offs and synergies) between multiple ES through time. Results Provision of ES Increased Through Time and Became More Variable Across Space. The provision of each service changed significantly

over time [test for the effect of time using space–time interaction analysis (STI), df = 393, P < 0.01]. The mean provision increased for almost every service, by between 20% and 94%, from 1971 to 2006 (Fig. 1 and SI Appendix, Fig. S1). Cultural services showed the greatest magnitude of change through time. Only the mean Significance Most approaches to quantifying and mapping ecosystem services (ES) focus on a single point in time. This static approach cannot provide insight into whether and how the provision of ES changes through time. We examined spatiotemporal ES dynamics by reconstructing the regional provision of nine ES over 35 y. Our approach demonstrated that individual services, ES bundles, and interactions among ES changed across both time and space. We also identified trajectories of ES bundle change and explained how these changes were driven by policy, biophysical, and socioeconomic characteristics. Our study demonstrates the limitations of assuming stationarity in ES and their relationships, and emphasizes the importance of taking into account both time and space in the assessment of multiple ES. Author contributions: D.R., J.M.R., and E.M.B. designed research; D.R. performed research; D.R. analyzed data; and D.R., J.M.R., and E.M.B. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. Data deposition: Ecosystem service estimates, socioeconomic and biophysical attributes for each municipality, for each date, are available on Dryad (dx.doi.org/10.5061/dryad.g4590). 1

To whom correspondence should be addressed. Email: [email protected].

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1502565112/-/DCSupplemental.

PNAS | October 27, 2015 | vol. 112 | no. 43 | 13411–13416

SUSTAINABILITY SCIENCE

Managing multiple ecosystem services (ES), including addressing trade-offs between services and preventing ecological surprises, is among the most pressing areas for sustainability research. These challenges require ES research to go beyond the currently common approach of snapshot studies limited to one or two services at a single point in time. We used a spatiotemporal approach to examine changes in nine ES and their relationships from 1971 to 2006 across 131 municipalities in a mixed-use landscape in Quebec, Canada. We show how an approach that incorporates time and space can improve our understanding of ES dynamics. We found an increase in the provision of most services through time; however, provision of ES was not uniformly enhanced at all locations. Instead, each municipality specialized in providing a bundle (set of positively correlated ES) dominated by just a few services. The trajectory of bundle formation was related to changes in agricultural policy and global trends; local biophysical and socioeconomic characteristics explained the bundles’ increasing spatial clustering. Relationships between services varied through time, with some provisioning and cultural services shifting from a trade-off or no relationship in 1971 to an apparent synergistic relationship by 2006. By implementing a spatiotemporal perspective on multiple services, we provide clear evidence of the dynamic nature of ES interactions and contribute to identifying processes and drivers behind these changing relationships. Our study raises questions about using snapshots of ES provision at a single point in time to build our understanding of ES relationships in complex and dynamic socialecological systems.

provision of cattle and flood control declined, by 20% and 30%, respectively. The STI analysis also showed that crop production, flood control, carbon storage and hunting activities were significantly different between municipalities at each time step (test for the effect of space using STI, df = 520, P < 0.01, Fig. 1). Additionally, the change in the provision of each service through time was not the same in all municipalities, as indicated by (i) the highly significant interactions between space and time for all services except flood control (STI, P < 0.01), and (ii) the increasing variation (i.e., SD) in the provision of each ES among municipalities.

1971 1981 1981 1991 2001 1991 2001 1971 STI (R2 = 0.04, p < 0.01*) S (R2 = 0.59 p < 0.01*) T (R2 = 0.09, p < 0.01*)

1971 1981 1991 2001 STI (R2 = 0.06, p < 0.01 *) S (p = 0.61) T (R2 = 0.13, p < 0.01*)

00

C stored/km2

11

22

Determinants of ES Dynamics Were Related to the Location of Municipalities. At all time steps, the spatial distribution of ES

-1

2 0

1

Pigs/km2

-1

-1

00

11

22

% Land under crop

Each ES Bundle Was Dominated by a Few ES. Cluster analysis partitioned the municipalities into seven groups based on the kind and amount of ES provided through time (Fig. 2). Bundles of services varied in the diversity of the ES provided (H, the effective number of ES, ranged from 4.35 to 7.14, Fig. 2). The

distribution across the landscape changed through time (Fig. 3). Although B2 (crops, cattle and flood control) was the dominant bundle in 1971 and 1976, provided by 41% and 44% of municipalities respectively, bundle types were more evenly distributed by 2006 (Fig. 3 A and B). Despite this increase in evenness, Moran’s I values showed that municipalities providing the same bundle became increasingly spatially clustered over time (in 1971 I = 0.11, P = 0.02; in 2006 I = 0.24, P < 0.01, SI Appendix, Table S1), forming a landscape-scale trade-off among ES bundles that is clearly visible in 2006 (Fig. 3C). This landscape-scale dynamic was due to changes in the bundle of services each municipality provided over time, which primarily followed four different trajectories (Fig. 3B). Municipalities providing B2 (crops, cattle and flood control) in 1971 changed to B3 (crop production) or B4 (carbon storage and flood control) by 2006, reflecting trajectories toward crop production specialization or toward rewilding following agricultural abandonment. Sixty-four percent of municipalities providing the B4 bundle in 1971 continued along the rewilding process through forest succession to provide high amounts of carbon stored and game animals (B6) by 2006. Twenty-seven percent of B5 municipalities in 1971 also changed to become B6 municipalities by the end of our study. This change occurred primarily in municipalities that were in the direct vicinity of the growing B6 cluster. The municipalities providing mainly recreational activities (B7) and animal production (B1) in 1971 continued to provide the same bundle of services through time. The number of municipalities providing the B1 bundle of services increased through time and expanded in spatial distribution into adjacent municipalities, resulting in the formation of a large cluster of municipalities that provided primarily animal production.

2001 1971 1981 1981 1991 1991 1971 2001 STI (R2 = 0.01, p < 0.01*) S (R2 = 0.55, p < 0.01*) T (R2 = 0.10, p < 0.01*)

2

All ES quantified rely on local ecosystems for their supply, including farm animals that are mostly fed with locally produced feed and outdoor recreation centers that involve outdoor activities in natural environments. Indicators were chosen based on the availability of primary data that was measured consistently through time and space. More detailed descriptions of the ES and historical data sources are provided in the SI Appendix. *Includes chickens, ducks, turkeys, and other poultry. † Flood control was quantified using the difference (i.e., δ) between the maximum and the mean number of flooding events over 5 y. ‡ Outdoor recreation centers are starting points for snowshoeing and skiing in winter or hiking trails in summer.

The Bundle Provided by Any Given Municipality Changed Through Time. The most common types of ES bundles and their spatial

1

Inhabitants/km Classes (0-7) km

2

0

Animals hunted/km2 Campsites/km2 Outdoor recreation centers‡/km2

1991 2001 1971 1971 1981 1981 1991 2001 STI (p = 0.34) S (R2 = 0.18, p < 0.01*) T (R2 = 0.18, p < 0.01*)

-1

δ(max − mean)†/5 y kg C/km2

Recreation centers/km2

Poultry*/km2 Cattle/km2 Pigs/km2 %

1.4

Provisioning Poultry Cattle Pork Crops Regulating Flood control Carbon sequestration Cultural Hunting Campsites Outdoor activities Covariables Population density Soil capability for agriculture Distance from Montreal

1

Indicators/units

0.6

Ecosystem services

bundle with the highest diversity value (bundle type 1, B1, H = 7.14) was characterized by multiple provisioning services with lower delivery of regulating and cultural services. The bundle with the lowest diversity value (B2, H = 4.35) had high delivery of provisioning services, but was dominated by crop and cattle production, with a high value for flood control. Bundles with intermediate diversity values typically had a few dominant ES and intermediate to low provision of other services. Among them, B3 was dominated by crop production and had intermediate values for cattle production and flood control. Bundles B4 and B5 were characterized by municipalities with intermediate production of crop and cattle but with other high-producing services (e.g., regulating services in B4, campsites in B5). Bundle type 6 (B6) characterized municipalities that provided high amounts of carbon storage and game animals, whereas municipalities in B7 were highly dominated by recreational activities.

Amplitude in the flood frequency

Table 1. Ecosystem services quantified from 1971 to 2006 across the Montérégie in Quebec, Canada, and covariables used to explain the spatial variability in the distribution of ES

1971 1971 1981 1981 1991 2001 1991 2001 STI (R2 = 0.13, p < 0.01*) S (p = 0.7) T (R2 = 0.31, p < 0.01*)

Fig. 1. Change in the provision of five selected ES through time (mean ± SD across all municipalities). Each ES was standardized to unit variance to allow comparison among the values. Results of the space–time interaction (STI) analysis performed for each ES are presented below the corresponding plot. The independent tests for spatial and temporal structures were abbreviated S and T respectively. We used α = 0.05, and the P value was calculated after 999 permutations. A significant STI indicates that the temporal change of a given ES was not the same in all municipalities.

13412 | www.pnas.org/cgi/doi/10.1073/pnas.1502565112

Renard et al.

relationships between animal production and cultural services shifted from trade-off, or no correlation in 1971 and 1976, to positive, significant, relationships by 2006. Relationships involving flood control were highly variable through time, showing no clear pattern of correlation with other services. Discussion We provide empirical evidence that the provision of individual ES and the relationships between ES are dynamic through time and space. Our analyses further showed that trajectories of change through time were not uniform across the region, but were instead related to the spatial distribution of environmental, social, and economic characteristics. We also showed that the relationships between services can shift through time. These results indicate that the commonly-held assumption that the provision of ES and their relationships stay the same over time is likely to be incorrect, and show the value of additional information that comes with understanding how complex social-ecological systems change through time. The provision of all but three of the ES that we quantified increased on average from 1971 to 2006. This trend is influenced by the types of ES we quantified in our study. In particular, we

Fig. 2. Ecosystem service bundles and the effective number of ES (H) provided in each bundle. Each petal in the bundles is associated with a symbol corresponding to an ES listed in Table 1. To facilitate comparison among services, ES abundances were normalized by the maximum ES value obtained for each bundle. The length of each petal is proportional to the relative abundances of the other ES within each bundle (petals are comparable within bundles).

ES Relationships Changed Over Time. At the regional scale, relationships among ES also changed through time (Table 2), in terms of both the type of relationship and its strength (strong: −0.5 ≤ r ≥ 0.5; moderate: −0.3 ≤ r ≥ 0.3; weak: −0.2 ≤ r ≥ 0.2). All cultural services showed a synergistic relationship with carbon storage but the strength of this relationship varied through time. While the relationship between hunting activity and carbon storage, already significant in 1971, became stronger through time (from 0.39 in 1971 to 0.73 in 2006, P < 0.01), the relationship between recreational activities and carbon storage became weaker (from a significant correlation coefficient of 0.24 in 1971 to a nonsignificant coefficient of 0.10 in 2006). Crop production showed consistent and mainly significant trade-offs with every cultural service and with carbon storage, whereas the Renard et al.

SUSTAINABILITY SCIENCE

bundles across the landscape was related to biophysical and socioeconomic attributes of the region (redundancy analysis applied over all time steps, R2 = 0.22, SI Appendix, Fig. S2). Municipalities providing primarily farm animals (B1), and carbon storage and game animals (B6) were distant from Montreal, had low population density, and soils with the lowest potential for crop production. Municipalities that specialized in recreational activities (B7) were located in areas with high population density, where soils also have low agricultural capacity. Areas where crops were produced (B2 and B3) were located where soils have the best potential for cultivation, at intermediate distances from Montreal, and with intermediate population densities. Municipalities providing B4 (carbon storage and flood control) and B5 (mix of food production and campsites) showed less clear relationships with explanatory variables. This finding aligns with our other results that showed that municipalities providing B4 and B5 changed through time (mostly replaced by B6, Fig. 3B), along with a change in their relationships with the explanatory variables.

Fig. 3. ES bundle dynamics over time (A and B) and across space (C). (A) The table shows the number of municipalities providing each bundle, for each time step from 1971 to 2006. (B) The web diagram shows the main trajectories of change that municipalities followed from one bundle to another between 1971 and 2006. The thickness of the lines is proportional to the number of municipalities in 1971, which follow the trajectory. (C) The maps display the spatial distribution of each bundle at each time step.

PNAS | October 27, 2015 | vol. 112 | no. 43 | 13413

Table 2. Pairwise correlations between ES through time

White indicates positive correlations that we defined as synergies (i.e., r > 0.1); dark gray indicates negative correlations that we defined as trade-offs (i.e., r < −0.1) and the star subscript indicates a significant relationship (P < 0.05); light gray reflects a relationship with no correlation (−0.1 < r > 0.1), always nonsignificant.

found very little data suitable for reconstructing regulating ES, which has almost certainly limited our ability to measure important environmental costs that come with increasing provisioning services (25). Our results also showed that the set of ES provided diversified, with cultural services, represented mainly by recreational services in our dataset, being increasingly provided across the region. However, the increasing provision of diverse services at the regional scale did not mean that each municipality provided a diverse bundle of ES. Instead, most municipalities specialized in a small set of services, resulting in seven unique bundles differentiating through time. Each bundle tended to be fairly specialized in one or a few services. Even the bundle with the highest diversity (B1) was primarily composed of provisioning services of various types, with low amounts of regulating and cultural services. The specialized nature of the bundles provided by municipalities in our study region, historically and today, is consistent with most bundles of ES identified in similar and other mixed-used landscapes (23, 24, 26). Although the provision of multifunctional bundles has proven to be possible at broader scales (26, 27), and is both socially desirable (28) and a target for enhancing ecosystem stability and human well-being (29), our results suggest that achieving multifunctional bundles, at least at small scales, has been difficult in the past. One reason could be that the diversification of bundle types, along with their specialization at the municipality scale, has produced economic benefits for large and specialized agricultural units (30). Reconstructing the trajectory of bundle changes through time can contribute to understanding how and why the provision of ES follows different trajectories, and identifying factors which might encourage greater multifunctionality in human-dominated systems. 13414 | www.pnas.org/cgi/doi/10.1073/pnas.1502565112

Our analyses also revealed that the different types of ES bundles have been increasingly structured in space in large clusters of contiguous municipalities. This spatial clustering has led to the recent emergence of a landscape-scale trade-off between provisioning (B1 and B3), regulating (B4 and B6), and cultural services (B7). Based on our analyses, we suggest that two underlying processes are responsible for the formation of such a trade-off: (i) the bundles provided through time by any given municipality tend to become more specialized; and (ii) the location of municipalities dominated this trend, leading to the provision of particular ES where biophysical conditions are most suitable or where people desired these services. Changes in the bundle of services provided by municipalities might also be explained by changes in the agricultural production context at provincial and international scales. Encouraged by an increase in the market value of corn, combined with the limited prospects for development in some agricultural sectors, Quebec adopted its first grain self-sufficiency policy in 1972 (31). The operationalization of this policy, supported by subsidies and advances in technology, encouraged the production of cash crops at the expense of the dairy industry and hayfields that had been dominant since the 1850s in our study region (31). The effects of these provincial policy changes were well reflected in two contrasting trajectories of ES bundles we identified between 1971 and 2006, particularly in the agricultural specialization (change from B2 to B3) and the rewilding after field abandonment (change from B2 to B4, and from B4 to B6 with forest succession). In parallel with provincial changes, a global trend toward intensification of pork production emerged and strengthened during the 1970s (32), leading to the spatial expansion of the bundle specialized in farm animal production (B1). A similar trend toward specialization of agricultural production has been described Renard et al.

Renard et al.

particular with provisioning services. Our ability to further understand how and why the supply of multiple ES and their relationships change through time can be much improved by promoting and designing long-term monitoring programs now for future ES historical work. Conclusions Our study contributes to a growing line of evidence showing that ignoring time in ecology can be perilous (6, 43–45). Just as today’s ecological communities are not considered a reliable representation of what communities may exist under future conditions (46), our study indicates that we should not expect ES supply to remain constant in time or space. Indeed, our spatiotemporal approach to analyzing multiple ES showed clear evidence of the dynamic nature of ES, their delivery, and their interactions. We have demonstrated that temporal dynamics in ES deserve more attention and hold potential to improve models of ES dynamics. Methods Quantification of ES. We quantified nine ES, including provisioning (n = 4), regulating (n = 2), and cultural services (n = 3). Table 1 describes the indicators used to estimate each ES and the sources of the data. Data Accessibility. Ecosystem service estimates, socioeconomic and biophysical attributes for each municipality, for each date, are available on Dryad (dx.doi. org/10.5061/dryad.g4590). Data used to quantify hunting activities, aboveground carbon storage, and part of pork production are protected by license agreements. See SI Appendix for data request procedure. Temporal and Spatial Scales. Each ES was quantified over 35 y, from 1971 to 2006, at 5-y intervals. This time period covered major policy changes in the study region (30). We quantified ES through time at the scale of administrative municipalities (n = 176) in the region of Montérégie in southern Quebec, Canada. Municipality area ranged from 27 to 256 km2 with an average of 80 km2. This spatial resolution is not only often used in historical records but is also a scale at which land-use management decisions are made. Because the boundaries of municipalities changed through time and the temporal resolution of historical data varied, we standardized the data to the same temporal and spatial scale. Details on the methodology used are presented in the SI Appendix. Data Set Analyzed. The final data set comprised 131 municipalities (44 municipalities were excluded either because there were no data recorded or because data were missing for at least four time steps), eight time steps, and nine ES. Missing data values were filled in with the average value of the given ES, over all of the other time steps, for the municipality of interest. This method has the advantage of limiting the weight of missing data values in multivariate statistical analysis but can obscure temporal trends. To test this, we compared the results we obtained to the results derived by filling in missing values using an interpolation method. We obtained very similar results. More advanced methods to handle missing data have been developed (47) but cannot yet be fully used for multivariate analysis. Before analysis, the final data set was transformed using the x’ = sqrt(sqrt(x)) transformation to meet assumptions of normality and standardized to unit variance and zero-mean to cope with the diverse units and ranges of variation of our data. All statistical analyses were performed using the software R v3.0.2 (48). Spatiotemporal Changes in the Provision of Individual ES. For each year, we calculated the mean and the SD of each of the nine ES, across all municipalities, to examine respectively the temporal trends in the provision of individual ES at the regional scale and the spatial variation of their provision through time. To statistically examine the STI, as well as the separate temporal and spatial changes in the provision of ES, we used the STI analysis following Legendre et al. (49). This method accounts for the absence of site replication in space–time sampling design. The time and space variables are coded using principal coordinates of neighbor matrices (PCNM) in a two-way analysis of variance (ANOVA) model 5 (999 permutations). A significant STI indicates that the spatial distribution of ES has changed through time or that the temporal trends were not the same in all municipalities. When STI was significant, we tested for separate temporal variation site to site and for spatial distribution time to time using a one-factor, ANOVA nested model. When STI was nonsignificant, we tested the significance of the space and time effects without replications using an ANOVA Model II. We applied the STI analysis to each ES individually.

PNAS | October 27, 2015 | vol. 112 | no. 43 | 13415

SUSTAINABILITY SCIENCE

in the United-States (33), in Europe (34), and elsewhere (35). Such specialization has often been paired with the development of capital- and technology-intensive agriculture, spurred by the economic efficiencies of operating at large scales and increased demands from international markets (36). Local environmental and social characteristics helped determine where changes in ES bundles happened across the region. For example, specialization in crop production happened primarily in municipalities with the best agricultural soils, whereas regulating services became predominant (B4 and B6) where soils were of limited use for crop production. Whereas farm animals were kept far from urban centers, the time series of maps of ES bundles through time showed that bundles comprising recreational activities and camping (B5 and B7) got closer to the main urban center through time. Mapping the current spatial distribution and congruence of ES has received much attention among the community of ES researchers over the last few years. These maps are used to target locations for management or conservation (37–39) and also to predict the future patterns of variation in the delivery of ES (40, 41). However, our results showed that such maps change through time, and that political decisions and local factors interact to shape these changes. Our study thus raises questions about the ability of an approach focused on mapping ES provision at one point in time to adequately support landscape management decisions aimed at long-term goals, which are by necessity rooted in the temporal dynamics of complex systems. Our analysis of temporal changes in pairwise interactions among ES revealed that interactions are not fixed in time. In our study region, animal production and cultural services (hunting activities in particular) shifted from a trade-off or no interaction at the start of the study period to an apparent synergistic relationship by the end. This change may be explained by changes in animal production methods from traditional breeding in extensive outdoor enclosures to production of animals inside specialized buildings (42). Although changes in farm management might have contributed to enhancing spatial compatibility with game animal movement and hunting activities, this likely led to more intense trade-offs with other services, such as water quality, erosion control and aesthetic appreciation for which historical data were not available. Relationships involving recreational activities also changed through time. The strength of the synergistic relationship between recreation and carbon storage decreased, whereas the strength of the trade-off between recreation and crop production decreased. These results mirrored changes in the spatial distribution of bundles comprising recreational activities and campsites (B5 and B7). At the beginning of the time period, these bundles were provided by the most remote municipalities, which were also the most forested. Through time, these bundles of services were increasingly provided closer to urban centers, in the direct vicinity of crop production areas. Although less associated with forests in these areas, activities associated with outdoor recreation centers, such as hiking, snowshoeing, or crosscountry skiing, may have been developed in wetlands and riversides. Changes in human demand for more accessible recreational services, or willingness and interest in recreating in agricultural areas, could explain the relaxation of the trade-off with crop production. With changing trade-offs, the possibility of new synergies can emerge between provisioning and cultural services, notably through agrotourism activities. Our study brings together the domains of ES science and historical ecology. We included ES that depend on natural ecosystems for their supply and were relevant to our study region. However, availability and quality of primary data, a challenge for all historical work, exerted a strong filter on the set of ES we could study and the indicators we could use. Ecosystem services for which long-term primary data are often unavailable (for example, pollination, water purification, aesthetic value, or sense of place) are usually not represented in historical studies like ours, leading to a reduced recognition of key trade-offs, in

Identification and Dynamics of ES Bundles. We identified ES bundles, i.e., mix of positively correlated ES provided together in the same place and at the same time, using a K-means clustering analysis on the entire time series. We selected the best partition based on the “simple structure index”. We analyzed the diversity of the set of ES provided in each bundle type using a transformation (H) of the Gini–Simpson’s entropy index (S): H = 1/(1 − S) (50). This transformation provides a measure of the “effective number of services,” which is the number of equally common services required to obtain the Gini–Simpson’s index. Using this transformation, diversity measurements have the same units and the same properties, no matter what was the diversity or entropy index used originally thus facilitating cross-study comparisons. We mapped the ES bundles for each year using ArcGIS (51) to visualize changes in their spatial distribution through time and determined their spatial clustering using Moran’s I (52). Finally, we calculated the percentage of municipalities changing from one bundle to another through time to examine the main temporal trajectories of change that municipalities followed between 1971 and 2006.

distance from the main urban center calculated based on digital boundary shapefile for 2006) and biophysical attributes (i.e., agricultural land capability, sources of data are detailed in the SI Appendix) of the region using a redundancy analysis (RDA). We controlled the RDA for temporal variation to specifically assess the link between the spatial distribution of the provision of ES and the characteristics of the region that did not change through time. The relationship was tested using a permutation test (53). Change in the Relationships Among Multiple Ecosystem Services Through Time. We performed Spearman correlations among each pair of ES (n = 26 pairs) for each time step to assess changes in the relationships among services through time.

Drivers of ES Dynamics Through Time and Space. We analyzed the relationship between the provision of ES and socioeconomic (i.e., population density,

ACKNOWLEDGMENTS. We thank C. Albert (University Aix-Marseille, UMR CNRS 7263/IRD 237), D. B. McKey (University Montpellier II, UMR CNRS 5175, Institut Universitaire de France), and two anonymous reviewers for their helpful comments; and Emily Clark (McGill University) for her language edits. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) in the form of a Strategic Projects grant (to E.M.B. and J.M.R.) and Discovery Grants (to E.M.B. and J.M.R.), as well as funds provided by the Trottier Institute for Science and Public Policy (TISPP).

1. Rodríguez JP, et al. (2006) Trade-offs across space, time, and ecosystem services. Ecol Soc 11(1):28. 2. Bennett EM, Peterson GD, Gordon LJ (2009) Understanding relationships among multiple ecosystem services. Ecol Lett 12(12):1394–1404. 3. Carpenter SR, et al. (2009) Science for managing ecosystem services: Beyond the Millennium Ecosystem Assessment. Proc Natl Acad Sci USA 106(5):1305–1312. 4. Nicholson E, et al. (2009) Priority research areas for ecosystem services in a changing world. J Appl Ecol 46(6):1139–1144. 5. Swetnam TW, Allen CD, Betancourt JL (1999) Applied historical ecology: Using the past to manage for the future. Ecol Appl 9(4):1189–1206. 6. Adler PB, Lauenroth WK (2003) The power of time: Spatiotemporal scaling of species diversity. Ecol Lett 6(8):749–756. 7. Foster D, et al. (2003) The importance of land-use legacies to ecology and conservation. BioScience 53(1):77–88. 8. Rhemtulla JM, Mladenoff DJ (2007) Why history matters in landscape ecology. Landscape Ecol 22(1):1–3. 9. Adler PB, et al. (2005) Evidence for a general species-time-area relationship. Ecology 86(8):2032–2039. 10. McKey D, et al. (2010) Pre-Columbian agricultural landscapes, ecosystem engineers, and self-organized patchiness in Amazonia. Proc Natl Acad Sci USA 107(17): 7823–7828. 11. Dupouey JL, Dambrine E, Laffite JD, Moares C (2002) Irreversible impact of past land use on forest soils and biodiversity. Ecology 83(11):2978–2984. 12. Fraterrigo JM, Turner MG, Pearson SM, Dixon P (2005) Effect of past land use on spatial heterogeneity of soil nutrients in southern appalachian forests. Ecol Monogr 75(2):215–230. 13. Hermy M, Verheyen K (2007) Legacies of the past in the present-day forest biodiversity: A review of past land-use effects on forest plant species composition and diversity. Ecol Res 22:361–371. 14. Francis R, Hare S (1994) Decadal scale regime shifts in the large marine ecosystems of the North-east pacific: A case for historical science. Fish Oceanogr 3:279–291. 15. Biggs R, Carpenter SR, Brock WA (2009) Turning back from the brink: Detecting an impending regime shift in time to avert it. Proc Natl Acad Sci USA 106(3):826–831. 16. Andersen T, Carstensen J, Hernández-García E, Duarte CM (2009) Ecological thresholds and regime shifts: Approaches to identification. Trends Ecol Evol 24(1):49–57. 17. Egan D, Howell EA (2001) The Historical Ecology Handbook: A Restorationist’s Guide to Reference Ecosystems (Island Press, Washington, DC). 18. Jackson ST, Hobbs RJ (2009) Ecological restoration in the light of ecological history. Science 325(5940):567–569. 19. Szabó P, Hédl R (2011) Advancing the integration of history and ecology for conservation. Conserv Biol 25(4):680–687. 20. Jiang M, Bullock JM, Hooftman DA (2013) Mapping ecosystem service and biodiversity changes over 70 years in a rural English county. J Appl Ecol 50(4):841–850. 21. MacDonald GK, Bennett EM (2009) Phosphorus accumulation in Saint Lawrence river watershed soils: A century-long perspective. Ecosystems 12(4):621–635. 22. Lautenbach S, Seppelt R, Liebscher J, Dormann CF (2012) Spatial and temporal trends of global pollination benefit. PLoS One 7(4):e35954. 23. Haines-Young R, Potschin M, Kienast F (2012) Indicators of ecosystem service potential at European scales: Mapping marginal changes and trade-offs. Ecol Indic 21: 39–53. 24. Raudsepp-Hearne C, Peterson GD, Bennett EM (2010) Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proc Natl Acad Sci USA 107(11):5242–5247. 25. Millennium Ecosystem Assessment (2005) Ecosystems and Human Well-Being: Synthesis (Island Press, Washington, DC). 26. Queiroz C, et al. (2015) Mapping bundles of ecosystem services reveals distinct types of multifunctionality within a Swedish landscape. Ambio 44(Suppl 1):89–101. 27. Foley JA, et al. (2005) Global consequences of land use. Science 309(5734):570–574. 28. Martín-López B, et al. (2012) Uncovering ecosystem service bundles through social preferences. PLoS One 7(6):e38970.

29. Pan Y, Xu Z, Wu J (2013) Spatial differences of the supply of multiple ecosystem services and the environmental and land use factors affecting them. Ecosyst Serv 5:4–10. 30. Chavas JP (2001) Structural change in agricultu ral production: Production economics, technology, and policy. Handbook of Agricultural Economics, eds Gardner B, Rausser J (Elsevier, North Holland, Amsterdam). 31. Domon G, Bouchard A, Gariépy M (1993) The dynamics of the forest landscape of Haut-Saint-Laurent (Quebec, Canada): Interactions between biophysical factors, perceptions and policy. Landsc Urban Plan 25(1-2):75–83. 32. Ruiz J, Domon G (2009) Analysis of landscape pattern change trajectories within areas of intensive agricultural use: Case study in a watershed of southern Québec, Canada. Landscape Ecol 24(3):419–432. 33. Hendrickson J, Sassenrath GF, Archer D, Hanson J, Halloran J (2008) Interactions in integrated US agricultural systems: The past, present and future. Renew Agric Food Syst 23(4):314–324. 34. Stoate C, et al. (2009) Ecological impacts of early 21st century agricultural change in Europe–a review. J Environ Manage 91(1):22–46. 35. Matson PA, Parton WJ, Power AG, Swift MJ (1997) Agricultural intensification and ecosystem properties. Science 277(5325):504–509. 36. Bowman MS, Zilberman D (2013) Economic factors affecting diversified farming systems. Ecol Soc 18(1):33. 37. Chan KMA, Shaw MR, Cameron DR, Underwood EC, Daily GC (2006) Conservation planning for ecosystem services. PLoS Biol 4(11):e379. 38. Turner WR, et al. (2007) Global conservation of biodiversity and ecosystem services. BioScience 57(10):868–873. 39. Egoh BN, Reyers B, Rouget M, Richardson DM (2011) Identifying priority areas for ecosystem service management in South African grasslands. J Environ Manage 92(6): 1642–1650. 40. Costanza R, Wilson A, Troy A, Voinov A, Liu S (2006) The Value of New Jersey’s Ecosystem Services and Natural Capital (New Jersey Department of Environmental Protection, Trenton, NJ). 41. Nelson E, et al. (2009) Modeling multiple ecosystem services, biodiversity conservation, commodity production, and tradeoffs at landscape scales. Front Ecol Environ 7(1):4–11. 42. Fraser D (2005) Animal Welfare and the Intensification of Animal Production. An Alternative Interpretation (Food and Agriculture Organization of the United Nations, Rome). 43. Cumming GS, et al. (2005) Are existing global scenarios consistent with ecological feedbacks? Ecosystems 8(2):143–152. 44. Raudsepp-Hearne C, et al. (2010) Untangling the environmentalist’s paradox: Why is human well-being increasing as ecosystem services degrade? BioScience 60(8): 576–589. 45. Tilman D, May R, Lehman CL, Nowak MA (1994) Habitat destruction and the extinction debt. Nature 371:65–66. 46. Williams JW, Jackson ST (2007) Novel climates, no-analog communities, and ecological surprises. Front Ecol Environ 5(9):475–482. 47. Honaker J, King G (2010) What to do about missing values in time-series cross-section data. Am J Pol Sci 54(2):561–581. 48. R Development Core (2008) R: A language and environment for statistical computing. 49. Legendre P, De Cáceres M, Borcard D (2010) Community surveys through space and time: testing the space-time interaction in the absence of replication. Ecology 91(1): 262–272. 50. Jost L (2006) Entropy and diversity. Oikos 113(2):363–375. 51. ESRI (2011) ArcGIS Desktop: Release 10. 52. Moran PA (1950) Notes on continuous stochastic phenomena. Biometrika 37(1-2): 17–23. 53. Legendre P, Legendre L (2012) Numerical Ecology, 3rd English Edition (Elsevier Science BV, Amsterdam).

13416 | www.pnas.org/cgi/doi/10.1073/pnas.1502565112

Renard et al.