Landscape complexity is associated with crop yields across a large temperate grassland region

Agriculture, Ecosystems and Environment 290 (2020) 106724

Contents lists available at ScienceDirect

Agriculture, Ecosystems and Environment journal homepage: www.elsevier.com/locate/agee

Landscape complexity is associated with crop yields across a large temperate grassland region

T

Paul Galperna,*, Jess Vickrucka,b, James H. Devriesc, Michael P. Gavina a

Department of Biological Sciences, University of Calgary, 2500 University Drive NW, Calgary, AB, T2N 1N4, Canada Agriculture and Agri-Food Canada, 850 Lincoln Road., Fredericton, NB, E3B 47Z, Canada c Ducks Unlimited Canada, P.O. Box 1160, Stonewall, MB, R0C 2Z0, Canada b

ARTICLE INFO

ABSTRACT

Keywords: Crop yield Ecosystem services Landscape complexity Semi-natural areas Perennial vegetation Canadian prairies

Establishing semi-natural areas within annual croplands can provide habitat for beneficial organisms and ecosystem services to crops through a spillover effect. However, this approach to increasing landscape complexity may have no effect on crops, or it may promote pests, weeds and other disservices that reduce productivity. An argument for changing landscape complexity may be more persuasive if it is associated with higher crop yields. Here, we examine regions that vary in their landscape complexity and, therefore, may also naturally differ in the potential for ecosystem services, disservices and crop yields. Specifically, we examine crop-growing districts in the Canadian province of Alberta to test whether the presence of more non-crop land covers has increased crop yields. Our dataset covered about one-quarter of the seeded area in Canada between 2012 and 2017 consisting of 10,069 records representing average field-level yields reported to a crop insurance provider. In total, we analyzed summary data for 250,000 km2 of seeded area for seven grain crops. Using a functional regression approach, we found evidence for a plausibly positive association between yield and the non-crop land covers found within and near fields in four of seven crops. Landscape complexity, therefore, represented a measurable yield benefit for farmers, although the variance in yield explained by the landscape was small. These findings suggest there may be a low risk of disservices to crops from non-crop land covers in this region. Our study adds support at a broad geographic extent for initiatives that restore perennial and other semi-natural vegetation in annual cropping systems and suggests that, in this temperate grassland region, their promotion (e.g., as carbon stores or as biodiversity refugia) may have no adverse effects for crop production.

1. Introduction Meeting the food needs of a growing global population will place significant demands on agriculture (Godfray et al., 2010), requiring the continued expansion and intensification of croplands (Egli et al., 2018; Laurance et al., 2014), while presenting risks for the conservation of biodiversity (Dalu et al., 2017; Karp et al., 2018; Kleijn et al., 2009; Titeux et al., 2016; Tubiello et al., 2015). Farmers may be able to minimize impacts on natural systems and improve crop productivity by encouraging ecosystem service provision in agricultural landscapes (Bommarco et al., 2013). One approach is to retain, restore or establish new semi-natural areas within the boundaries of a field or in the neighbourhood of a crop (Fahrig et al., 2011; Landis, 2017). More complex agricultural landscapes (i.e., with more semi-natural or non-crop land covers) can host more abundant and diverse communities of organisms that provide ecosystem services to agriculture.

⁎

These beneficial organisms include flower-visiting insects, birds and the natural enemies of crop pests that use the greater number of habitats that may be available, leading to “spillover” into nearby crops (Chaplin‐Kramer et al., 2011; Duarte et al., 2018; Ekström and Ekbom, 2011; Gallé et al., 2019; Landis et al., 2000; Letourneau et al., 2011; Tscharntke et al., 2005; Vickruck et al., 2019; Winqvist et al., 2011). Increasing agricultural landscape complexity may go further than creating habitat for beneficial organisms. For example, encouraging perennial vegetation can provide carbon storage opportunities (Hungate et al., 2017; Lamb et al., 2016; Smith, 2014; Williams et al., 2018), while also conserving habitat for other biodiversity (Fischer et al., 2014; Phalan et al., 2011; Ponisio et al., 2016; Tscharntke et al., 2012). While there is potential for a suite of ecosystem services to be promoted through increased landscape complexity, implementation will depend on understanding the economic costs and benefits of these

Corresponding author. E-mail address: [email protected] (P. Galpern).

https://doi.org/10.1016/j.agee.2019.106724 Received 22 June 2019; Received in revised form 3 October 2019; Accepted 21 October 2019 0167-8809/ © 2019 Elsevier B.V. All rights reserved.

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

interventions (Benton et al., 2003; Bommarco et al., 2013). A persuasive metric for farmers may be crop yield. For example, the case for retaining an uncultivated field margin may be more compelling if crop yields can be shown to benefit from the ecosystem services that are hosted there. Equally, if uncultivated areas are found to have no negative impact on crop yields, farmers may be able to justify them as biodiversity refugia, or for carbon storage and other ecosystem service benefits. Research across a range of cropping systems have shown that diversifying on-field vegetation can improve crop yields and provide ecosystem services (Duarte et al., 2018; Garbach et al., 2017). Off-field, at the landscape scale, the potential for landscape complexity to influence crop yields has proven harder to measure (Bommarco et al., 2018; Garibaldi et al., 2011). Studies examining a variety of beneficial arthropods including spiders, carabid beetles, bees and other flower visitors, underline how proximity to natural features can increase pollination and pest control services as well as crop yields (Gallé et al., 2019; Garibaldi et al., 2011). Field size, which is likely to influence landscape complexity, habitat heterogeneity and the proximity of offfield habitats to crops (White and Roy, 2015) has been linked to reduced wind erosion of soils, pollination and improved crop yields in multiple pollen-limited crop systems (Garibaldi et al., 2016; Gene et al., 2019). Landscape complexity may also have a negative effect on crops. For example, semi-natural areas close to fields may be reservoirs for pests and weed seeds. Although meta-analyses suggest that insect pest abundances are generally unrelated to landscape complexity on average across a variety of cropping systems, the impact of these pests may be a consequence of complex population dynamics among the focal species, its host plants and its natural enemies, all of which can be influenced by landscape complexity (Chaplin‐Kramer et al., 2011; Rusch et al., 2016). Equally, weed plant distributions may be unrelated to landscape complexity, although this may vary with dispersal and other life history attributes of the plant (Alignier et al., 2012). Complex landscapes can also present risks for conservation objectives, for example, by providing pathways for invasive species to spread (O’Reilly-Nugent et al., 2016). While evidence suggests that ecosystem disservices to agriculture associated with landscape complexity should not be prevalent, in some systems and growing seasons, there remains potential for non-crop land covers to be a net drain on crop productivity. This body of work creates the expectation that crop yields can be improved where there is habitat within or near fields that can support diverse and abundant assemblages of organisms that provide services to agriculture, or where field sizes and land covers promote environmental conditions that support crop productivity. Data demonstrating this relationship have typically been obtained at the scale of one or several fields (but see Garibaldi et al., 2016). It remains unknown whether a positive or negative pattern may emerge at a larger spatial and temporal extent, across which the assemblages of species providing ecosystem services and disservices may vary, and the structure and composition of the habitats hosting these organisms is diverse. The aim of this study is to estimate the strength and the direction of the relationship between landscape complexity and yield across a large geographic extent. Following existing research, we hypothesize that the land cover diversity created by landscape complexity will influence the provision of both ecosystem services and disservices to crops (e.g. by leading to the spillover of organisms from nearby habitats, or by changing soil erosion potential, soil moisture and other local environmental conditions). This leads to the prediction of an association between crop yields and landscape complexity, with the direction of this association reflecting the net influence of services or disservices. We test this prediction by comparing crop yields in neighbouring administrative districts that differ naturally in their landscape and biological diversity, using measures of landscape complexity defined by heterogeneity in land use (i.e. whether the land is used for crops or otherwise). We estimate this relationship for seven of Canada’s major

grain crops over six consecutive years using a functional regression approach. 2. Materials and methods 2.1. Study system Our geographic focus is cropland in the province of Alberta, Canada, which represents about one-quarter of the seeded area in Canada (approx 104,000 km2 seeded; Statistics Canada, 2018). There is considerable diversity in land cover across Alberta’s agricultural zone. It includes the Mixed Grassland, Aspen Parkland, and Peace Lowland ecoregions of the Canadian prairies (Ecological Stratification Working Group, 1995) and is one of the most intensively cropped regions in the world (Foley et al., 2005). Dominant agricultural land uses include cropland (predominantly for cereal grain and oil-seed production), and introduced and native perennial grass forage lands (pasture and hay) for cattle ranching. Semi-natural areas within and adjacent to croplands that are not used for agriculture are dominated by perennial grasses and shrubs, with more tree cover at higher latitudes (Ecological Stratification Working Group, 1995). Much of this region is also characterized by the presence of numerous wetlands, many of which are small “pothole” basins embedded within crop fields (Doherty et al., 2018). Thus, the study region provides a broad gradient in land cover over which to test for a relationship with crop yield. Yield data were obtained from a large multi-annual agricultural data set for Alberta, providing average yield information summarized by municipal district, year and crop variety for seven of the region’s most important grain crops. These data were collected by the Agriculture Financial Services Corporation (AFSC), a state-owned enterprise providing crop insurance to the majority of farmers in Alberta. Municipal estimates of yield provided by AFSC were based on reports of harvests in each insured field and a precise accounting of the seeded area (i.e. excluding non-crop land cover within field boundaries). 2.2. Landscape complexity 2.2.1. Land cover mapping A land cover map for the agricultural regions of Alberta, Canada at 30 m resolution was developed by combining data from multiple published data sources (Table A1). It was then simplified thematically to a binary map (i.e. in crop production or otherwise). This classifies all parts of the landscape where crops are not present as contributing equally, regardless of the class of land cover they represent (i.e., noncrop land covers). The resulting product can be understood as a map of land use for crops, which is more interpretable across a large geographic region where there is a broad range of natural and semi-natural vegetation classes. Non-crop land covers included natural features that are common within and near Alberta fields, however the thematic resolution of the available data products did not permit their separation into consistent sub-classes for this analysis. These features included pasture, pothole wetlands, tree and shrub patches, uncultivated grassland patches, shelterbelts, and road-side grass margins all of which may influence ecosystem services and disservices to nearby crops (Asbjornsen et al., 2014). Paved and other human-modified areas such as roads, buildings, farm yards and in-field oil and gas well pads were included within the definition of non-crop land covers. These features may have some association with ecosystem services and disservices by introducing heterogeneity into the landscape; for example, by altering local environmental conditions, or because they may be surrounded by semi-natural land cover (e.g. weedy grasses, forbs and shrubs), and host interstitial micro-patches that act as habitat reservoirs for beneficial or pest organisms (Forman, 2009). The accuracy of the binary land cover map was assessed by randomly positioning ten 400 ha squares in each of four 3600 km2 scenes 2

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

multiple spatial scales. Unlike measuring landscape complexity using metrics such as class diversity, class evenness, shape, patch size or fragmentation (e.g., Boesing et al., 2017; McGarigal et al., 2012), this method does not require the specification of an analysis radius a priori, which may be beneficial when the distance at which non-crop land cover exerts an influence is unknown. The expected distribution for each municipal district was estimated as a continuous function of the distance from the centre of fields using a type of functional data analysis (FDA) known as function-on-scalar regression (Kokoszka and Reimherr, 2017; Wood, 2017). In these models, a continuous function Pi (s ) , is estimated from discrete data (e.g. Fig. 1, B, light line) giving the proportion of non-crop land cover at each annulus radius, s , in field, i . In order to find the mean function for fields in each municipal district, we fitted a model of the general form,

from across the province that were selected to provide a latitudinal gradient. All non-crop land cover objects visually identifiable on 0.5 m resolution true-colour snow-free satellite imagery (DigitalGlobe, 2016) were manually digitized and their areas compared to those on the land cover map. 2.2.2. Land cover within and near fields Fields across much of Alberta are regular in shape and size and typically form a portion of a square 259 ha (1 mile by 1 mile) section, a consequence of consistent survey practices in the nineteenth century (Larmour, 2005). Sections, or quarter sections within, are often planted in a single crop, and are surrounded by roads, tracks or other non-crop land covers that can be used to delineate field boundaries. The binary land cover map was processed using a clump detection algorithm (raster package for R; Hijmans, 2015) to find all such delineated crop patches, and the centroid and area of each patch was used to estimate the centre and size of a field respectively. Focal field centres for this analysis were randomly sampled from all candidate field centres with the constraint that field centres be a minimum of 3.5 km apart. This was done using an iterative algorithm that tested field centres for inclusion by drawing them from a randomly-permuted list, stopping when no more that met the constraint could be added. Land cover was measured at multiple distances from focal field centres by creating 52 annuli of increasing radii up to a maximum radius of 1514 m (illustrated in Fig. 1, A). The radii were chosen such that each annulus would have approximately the same total area to minimize bias associated with a smaller sample of land cover. This required that the distance between radii progressively decreased (Fig. 2, B, dark line). The proportion of non-crop land cover in each annulus was then measured from the land cover map ensuring cells that fell in more than one annulus were counted only in the smallest annulus in which they were measured.

Pi (s ) = µ (s ) +

mi

+

mi (s )

+

r 1 r 2 (r1i ,

r2i, s ) +

q (qi ,

s) +

i (s )

(1) where µ (s ) is the overall mean land cover function, mi is a constant effect and mi (s ) is the mean land cover function for the municipal district mi , and i (s ) is the residual function. The model assumed a Binomial datagenerating process and was fit using the refund package for R (Huang et al., 2016) which provides appropriate settings for FDA when interfacing with the mgcv generalized additive modelling (GAM) package for R (Wood, 2017). Non-linearities in the mean municipal district functions were estimated non-parametrically using penalized splines. Additional regressors included a three-dimensional tensor product for the term, r 1 r 2 (r1i, r2i, s ) to smooth spatial variation across municipal district eastings, r1, and northings, r2 . A two-dimensional tensor product for q (qi , s ) , represented the log -transformed area, q, of the field, and was used to control for differences among municipal districts in their field size distribution. The maximum number of knots that could be used for estimating municipal district penalized splines was set at 42 as a balance between overfitting and modelling abrupt changes in non-crop land cover that were evident across annulus distance classes (e.g. in the transition between field and field boundary). After fitting, estimated mean land cover functions were retrieved by setting the other covariates at each of the municipal district averages. These functions, each one representing expected land cover variation in a municipal district, were used as explanatory variables in the crop yield models.

2.2.3. Modelling regional landscape complexity We modelled landscape complexity by estimating the expected distribution of non-crop land covers continuously from field centres for each municipal district. This approach can be understood both as a metric for the proximity of non-crop land cover to a field centre, and more generally as an index of the variation in non-crop land cover at

Fig. 1. Illustration of the method used for the measurement of landscape complexity at multiple spatial scales. (A) A portion of the binary land cover map. (B) A curve describing landscape complexity (light line) was estimated by finding 52 annuli of increasing radii from the centroid of contiguous area of crop (a “field”) such that each annulus contained approximately the same number of cells (dark line). The number of non-crop cells was counted in each annulus (light line) and this process repeated for all fields sampled. Cells that fell in more than one annulus were counted in the smallest annulus in which they were measured.

3

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Fig. 2. Estimation of landscape complexity at multiple spatial scales using non-crop land cover measured continuously from field centroids. (A) Locations of focal fields (n = 8507; points) in municipal districts in Alberta, Canada (n = 66; regions). (B) Over-plotted semi-transparent Gaussian kernel density estimates of the distribution for municipal districts at three selected annuli show considerable variation both within and among districts and that this distribution varies with spatial scale. (C) Curves estimated from these data using function-on-scalar regression represent the expected proportion of non-crop land cover in a municipal district at multiple spatial scales. Each curve is used to summarize landscape complexity in a single district, and serves as a predictor in crop yield analyses. These have been coloured using the latitude of the district centroid.

2.2.4. Estimating regional field size We estimated the field size distribution in each municipal district to provide interpretation for the annulus radius and to demonstrate its relationship to conditions on the ground. We used the distance from focal field centres to the nearest road on our land cover product to estimate the smallest dimension of a field. Following the same procedure as Section 2.2.2, the number of road pixels in each annulus was determined and the smallest annulus with at least 3 cells of road (90 m length) was recorded as the field boundary. This threshold was applied to minimize misclassifying tracks within fields as boundaries. Distributions of the distances to the nearest road were visualized for each municipal district using the default Gaussian kernel density estimator settings provided by the ggplot2 package for R (Wickham, 2016).

used to model random effects and spatial covariation (Wood, 2017). These included a random intercept of crop variety, (vk ), to represent differences among varieties, vk , in their potential yields, and a Markov random field smooth function, (mi ) , which modelled covariation among municipal districts adjacent to district mi . This latter random field was intended to capture the effects of otherwise unmodelled, but spatially-variable, environmental conditions such as soil or climate, as these are expected to demonstrate geographic patterns at the extent of this analysis (Ecological Stratification Working Group, 1995) and therefore drive covariation among yields in spatially-adjacent municipal districts. Two models were fit for each crop assuming a Gaussian or a Scaled t process, as heavy tails were evident in initial model checks. The model of each pair with the lowest AIC was reported. To assess the predictive importance of the non-crop land cover term we also fitted models with this functional term excluded (reduced models). In addition to including smooth terms intended to control for spatial and temporal variation, we inspected model residuals, i , for evidence of autocorrelation. Moran correlograms at twenty distance lags with a permutation test (n = 1000) at each were used for spatial autocorrelation (ncf package for R; Bjørnstad, 2016). For serial autocorrelation, the estimated mean and variance of the distribution of one- and two-year autocorrelations for each municipal district-crop variety combination were compared to values from distributions of permuted data (n = 1000).

2.3. Estimation of crop yield associations 2.3.1. Modelling crop yield The association between crop yield and non-crop land cover was modelled at the municipal district level using a second class of FDA model known as scalar-on-function or signal regression (Kokoszka and Reimherr, 2017; Marx and Eilers, 1999; Wood, 2017). In these scalar response models, crop yield, yijk , in each municipal district, i , year, j , and crop variety, k , was modelled using the area under the mean land cover functions, Pi (s ) , retrieved from the function-on-scalar model (1). We fit models of the general form,

yijk =

(s ) Pi (s ) ds +

(aj ) +

(r 2 i ) +

(v k ) +

( mi ) +

ijk

2.3.2. Exploration of effect size We mapped the predictions of the yield model to illustrate the effect on crop yield of redistributing crop land cover from across a municipal district to non-crop land cover within the boundaries of a typical field. The fitted model was used to estimate the change in yield in each municipal district associated with modifying the shape of the land cover function. The mean function, Pi (s ) for each municipal district, i , was transformed by adding to the proportion of non-crop land cover at a subset of annulus radii such that the area added corresponded to a 1 % reduction in the total crop area predicted for each municipal district. Estimates for these maps were made while controlling for crop variety and year.

(2)

where ijk are residuals, and the coefficient function gives the estimated effect of land cover variation on yield, which is the primary effect of interest. Yield, yijk , was centred and scaled to permit comparisons of parameters among crops. Models were fit using restricted maximum likelihood estimation also with the refund package as an interface to mgcv (Wood, 2017). Additional covariates modelled the structure of the yield data. A smooth year term, (aj ) , modelled year, aj , as a penalized spline to control for the expected variation among years in crop productivity. A smooth northing term, (r 2i ) , as a penalized spline was added to test for a latitudinal productivity gradient. Splines were also 4

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Fig. 3. The distribution of distances between field centres and the nearest road by municipal district based on a Gaussian kernel density estimate of all sampled fields in each district. Expected modes, given field regularization by nineteenth century land surveys, are shown using dotted vertical lines.

3. Results

3.2. Landscape complexity

3.1. Agricultural data

Assessment of error on the land cover map revealed relatively small differences between manual object identification and the mapping product used in these analyses (maximum 4.4.% misclassification of non-crop as crop, and maximum 4.8.% misclassification of crop as noncrop; Table A3). Estimates of landscape complexity were based on the assessment of 8507 field centres in 66 municipal districts (Fig. 2, A). The distribution of non-crop land cover differed within and among districts (Fig. 2, B). Most districts had significant non-linearities in their estimated functions of land cover over distance, implying that the

The agricultural data set included yield data from approximately 250,000 km2 of crop fields between 2012 and 2017, and was aggregated by year, municipal district (n = 66), and crop variety (n = 420) for canola (a type of oilseed rape), barley, oat, pea and three classes of wheat (n = 10,069 records; Table A2).

5

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Fig. 4. Estimated coefficient functions and their 95 % confidence regions from scalar-on-function regressions of non-crop land cover at multiple spatial scales (i.e. landscape complexity functions) on yield (scalar) for seven crops. Coefficient functions must be interpreted in their entirety, and here imply that effects of non-crop land cover on yield represent a balance between the amount recorded at small and large spatial scales (i.e. annuli).

distribution of non-crop cover with respect to field centres had a discernable pattern across fields within that district (Table A4, a; i.e. that they differed in their landscape complexity). There was also considerable variation in the estimated mean amount of non-crop land cover in a district (Table A4, b). The predicted functions from this model demonstrate that the amount and distribution of non-crop land cover varies considerably across the agricultural region, and that latitudinal gradients (e.g., in moisture or productivity) are not highly correlated with this variation (Fig. 2, C).

that fields have been regularized by nineteenth century surveys on 1 mile by 1 mile sections (a grid spacing of approximately 1600 m by 1600 m). In many districts there was a mode at, or less than, 800 m, implying that annulus distances greater than this can generally be interpreted as “outside the focal field”. 3.4. Estimation of crop yield associations The proximity of non-crop land cover to the centre of a field was significantly associated with the mean yield per unit area in all seven crops (Fig. 4; Table 1). Examining coefficient functions, we found evidence that for barley and hard red spring wheat, at distances closer to the centre of the field, the proportion of non-crop land cover had a positive effect on yield. For canola, pea, oat, and Canadian prairie spring wheat, non-crop land cover had the most plausibly positive

3.3. Field size We found that the distribution of the distance between field centres and the nearest road had two or three modes in most municipal districts (400 m, 800 m and 1200 m; Fig. 3), corresponding to the expectation 6

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Table 1 Modelling of crop yield associations with non-crop land cover at the municipal district level using several smooth predictors that included a functional (signal-spline) predictor to represent non-crop land cover variation at multiple spatial scales (i.e. landscape complexity). Parameter estimates for seven crop models are presented.

CANOLA

BARLEY

PEA

OAT

WHEAT (Hard Red Spring)

WHEAT (Canadian Prairie Spring)

WHEAT (Soft White Spring)

Parameter

Parameter type

edf

P

σparam

ICC

Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.57) ΔAICfull-reduced = −1.74 Scaled t (v = 5.85)

spline signal-spline spline random effect Markov random field

2.72 2.88 4 92.21 48.84

0.241 < 0.001 < 0.001 < 0.001 < 0.001

0.04 7.26 0.35 0.12

0.003 0.994 0.279 0.043

1.01 2 3.99 34.04 45.32

0.981 < 0.001 < 0.001 < 0.001 < 0.001

< 0.01 6.63 0.47 0.16

< 0.001 0.991 0.361 0.058

3.24 2.72 3.97 8.98 25.94

0.0782 0.0131 < 0.001 < 0.001 < 0.001

0.04 6.7 0.17 0.76

0.002 0.987 0.048 0.021

1.22 3.09 3.94 6.22 20.74

0.0477 < 0.001 < 0.001 < 0.001 < 0.001

0.24 0.07 5.53 0.29 0.7

0.104 0.009 0.984 0.149 0.021

spline

2.37

0.509

signal-spline spline random effect Markov random field

2 3.99 37.11 48.2

< 0.001 < 0.001 < 0.001 < 0.001

< 0.01 7.27 0.36 0.15

< 0.001 0.995 0.309 0.075

spline

1

0.564

signal-spline spline random effect Markov random field

2.82 3.98 7.54 28.8

< 0.001 < 0.001 < 0.001 < 0.001

0.06 8.44 0.27 0.18

0.01 0.995 0.171 0.084

spline

2.13

0.00274

0.7

0.738

signal-spline spline random effect Markov random field

2 3.95 0.36 11.48

0.0205 < 0.001 0.196 < 0.001

< 0.01 7.42 0.05 0.13

< 0.001 0.997 0.014 0.091

Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.62) ΔAICfull-reduced = −2.8 Scaled t (v = 7.6) Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.76) ΔAICfull-reduced = 2.78 Gaussian Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.67) ΔAICfull-reduced = −5.64 Gaussian Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.54) ΔAICfull-reduced = 1.72 Scaled t (v = 9.68) Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.60) ΔAICfull-reduced = −5.26 Scaled t (v = 11.37) Northing f (Non-crop, Distance) Year Crop variety Municipal district Scale (σЄ = 0.42) ΔAICfull-reduced = 0.64 Scaled t (v = 3.79)

n = 4101

2 R adj = 0.52

spline signal-spline spline random effect Markov random field

n = 2111

2 R adj = 0.47

spline signal-spline spline random effect Markov random field

n = 719

2 R adj = 0.42

spline signal-spline spline random effect Markov random field

n = 464

2 R adj = 0.51

n = 2082

2 R adj = 0.63

n = 452

2 R adj = 0.56

n = 140

2 R adj = 0.65

7

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Fig. 5. Exploration of the effect size associated with the estimated coefficient functions for four crops. These predicted maps show the proportional increase in yield associated with redistributing 1 % of the crop land cover from across a municipal district to non-crop land cover within or near the boundaries of a typical field. These maps are intended only to illustrate the effect size and do not account for other more important sources of variation.

effects when these covers occurred 600−800 m from field centres. In all cases, additional non-crop land cover at greater distances tended to have plausibly negative effects on yield (e.g. intercept at ∼1100 m for canola; Fig. 4). Coefficient functions from scalar-on-function regressions should be interpreted jointly across all values, implying that trends across the modelled range of the function are collectively associated with the response variable (Kokoszka and Reimherr, 2017). In the fitted models, the existence of both positive and negative coefficients implies that the effect on yield represented a balance between the amount of non-crop land cover near and far from the centre of the field. This balance, and therefore the magnitude of the effect on yield, varied

among the crops tested. For most crops, plausibly positive effects occurred at annuli less than 800 m (Fig. 4), a distance that could be considered within or near the boundaries of the majority of fields in the study region (Fig. 3). The most plausibly negative effects on yield were found much further from field centres (> 1000 m). Based on field size modeling, this distance would often fall beyond the boundaries of a typical field and consequently describe the landscape context. All crops demonstrated significant year-to-year variation in yield (Fig. A1), and there was evidence of a random effect of crop variety, implying differences among varieties in their baseline yields (Table 1). There was no evidence of a latitudinal effect on yield in five of seven

8

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

crops (Northing; Table 1). However, Markov random field smoothers in all crops indicated that neighbouring municipal districts covaried, implying that spatially-varying environmental factors such as climate and soil conditions that are likely to produce yield differences among districts, but are not explicitly included as variables, have implicit representation in the model (Table 1). Fig. A2 maps the smoothed Markov random fields for each crop and can be understood to approximate the effect on yield of these environmental factors. We found no evidence of spatial or serial autocorrelation in the residuals of any crop model.

with positive and negative associations emerging at different scales (Fig. 4). These results underline that field-level studies showing relationships between landscape complexity and ecosystem service provision to crops (Duarte et al., 2018; e.g., Gallé et al., 2019) are detectable at an aggregated district-level of analysis, and that relationships hold across a large agricultural region over multiple years, crops and varieties. Our results also suggest that increasing landscape complexity across a district can make a measurable contribution to average yields, but it is likely to be small relative to the variety planted, and factors beyond the control of growers such as year-to-year variability and spatial variation in soil and climate. If we assume constancy in these variables, our models predicted that changing landscape complexity by adding small amounts of land covers within or near the boundaries of fields can improve yields in aggregate and potentially compensate growers for the loss of cultivation area (Fig. 5). However, the practical impact of this observation is low given the small amount of variance in yield represented by non-crop land cover (e.g., ICC 0.01 depending on the crop; Table 1). As such, we advocate that our results be cautiously interpreted as providing evidence that the effect of landscape complexity on yield is plausibly positive, and recommend that the effect be explored at the field level in order to determine how the characteristics of non-crop land cover (e.g., class, area, and shape) contribute to yield, and to estimate the limits of this effect.

3.4.1. Prediction of crop yield There was a large difference among modelled variables in their contribution to variation in crop yield. In all crops the proximity to noncrop land cover emerged as the smallest variance component ( param , Table 1) and had the lowest correlation with variation in yield (Intraclass correlation coefficient, ICC; Table 1). These measures imply that land cover relationships with crop yield represented a very small proportion of the expected variation in yields. Indeed, in three crops AICfull reduced > 0 , demonstrating that the full model with the non-crop land cover functional term included had lower predictive power than the same model with this term dropped (Table 1). Ranked in terms of decreasing ICC, Canadian prairie spring wheat followed by oat, canola and pea were crops with the most evidence for covariation between yield and land cover (Table 1). Inter-annual variation emerged as the largest correlate of the observed variation in crop yield, with crop variety and municipal district Markov random field smoothers typically ranked second and third depending on the crop ( param , ICC; Table 1).

4.1. Interpretation of models A key outcome from the modelling is that the yield per area in a municipal district increases when fields in that district tend to have more non-crop land cover within or near their boundaries (i.e. < 800 m; Fig. 4). This statement requires a caveat because only part of the coefficient function is positive. The positive association closer to the centre of the field is discounted by the amount of non-crop land cover outside the boundaries of those fields (i.e. > 1000 m; Fig. 4). For example, in districts where many fields may be surrounded by pasture lands or forest the model suggests that the high amounts of non-crop land cover in the landscape context would reduce any positive yield effect of non-crop patches within the fields. It is easier to explain the main result if we understand the distance of the annulus from the field centre as a measure of spatial scale, rather than as an indicator of distance. First, recall that the non-crop land cover functions (Fig. 2, C) used as predictors of district yields are averages over focal fields within a district. As such, they describe the expected pattern of non-crop land cover moving outwards from the centre of the average field. Higher distance classes (> 1000 m) could be described as capturing the landscape context because annuli will typically fall outside the boundaries of a focal field and will therefore summarize land cover from several neighbouring fields and natural landscapes. Where there is less non-crop land cover at this scale, this implies that croplands cover more of the region, and production could be said to be extensive. In contrast, distance classes within or near the boundaries of typical fields (< 800 m) summarize the scale at which management interventions can feasibly occur; e.g., where there is less non-crop land cover at this scale, production could be described as more intensive. We therefore interpret each function (Fig. 2, C) describing variation in non-crop land cover as a multivariate index of landscape complexity integrated over multiple spatial scales. The fitted coefficient functions in the yield models (Fig. 4) also retain the multiple spatial scale interpretations. From this perspective,

3.4.2. Exploration of effect size Fig. 5 maps the predicted effects for each municipal district of converting 1 % of the crop area near an average field to non-crop land cover (i.e. 600 m–750 m from the centre of that field), while controlling for all other sources of yield variation. In all municipal districts with yield data for that crop there was at least a 1.1.% increase in the yield measured on a seed weight per unit area basis ¯

¯

¯

¯

( canola = 2.0%, barley = 1.4%, Cps wheat = 4.6%, oat = 5.1%). These four crops were those where the inclusion of the proximity to non-crop land cover functional term improved prediction (i.e., AICfull reduced < 0 ), and the distances chosen for conversion were those where the predicted non-crop land cover functions were plausibly above zero (i.e., as determined from Fig. 4), permitting illustration of the largest positive effects on yield. 4. Discussion Our analyses have implications for increasing landscape complexity, and show that yields are positively associated with non-crop land covers, although other environmental factors beyond the control of farmers are likely to be much more important drivers of crop productivity. Further, they demonstrate that any ecosystem disservices to crops associated with non-crop land covers are unlikely to be dominant, reducing disincentives for their conservation as habitats for biodiversity and as carbon stores. We used existing gradients in the amount and distribution of noncrop land covers across Alberta, Canada to test this association. The main finding was that mean crop yields in a district were associated with the mean proximity of non-crop land to field centres (Table 1),

9

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

higher yielding municipal districts are those with a certain pattern of landscape complexity. Interpreted in these terms, the coefficient function shows that districts where the landscape context has extensive production, but fields are relatively less intensive, will have a higher yield. The pattern at the landscape context scale (> 1000 m) appears to contradict the prediction that non-crop land cover will increase yields (Fig. 4). This negative association can be explained by the confounding effect of other environmental factors. Farmers will have an incentive to put more land into crop if that land tends to give them higher yields; for example, the climate and soil conditions in a district with more extensive production may be closer to the optimum for these crops, resulting in higher profitability and therefore more investment in clearing of non-crop covers. From these economic drivers alone, non-crop land cover should be expected to negatively correlate with other environmental factors that increase crop yield, and therefore extensive production should be positively associated with crop yield. The value of functional modelling is that it has allowed us to control for this effect of extensive production in the landscape context and its potential to mask an association between non-crop land cover and yield. We therefore advocate that future analyses of crop yield and land cover at an aggregate level should include measures of extensive production. To help readers understand the magnitude of the statistical effects represented by these functional regression models, we examined the implications of changing the distribution of non-crop land cover in a district. Projections using this model show that removing 1 % of the cropland in a municipal district from production will result in at least a 1 % increase in yields measured on an areal basis, provided this removal occurs within or near the boundary of a typical field (i.e. 600–750 m; Fig. 5). While caution is certainly appropriate with such extrapolations, these figures illustrate the magnitude of the intensive production effect, and that it is potentially large enough to compensate for the losses in total crop area that would occur through increasing landscape complexity. While the estimated coefficient functions show there is potential for yield improvements, the variance associated with non-crop land cover was very low in all crop types, and in pea, hard red spring wheat and soft white spring wheat crops did not improve the model fit ( AICfull reduced > 0 ; Table 1). This underlines that at the aggregate scale of this analysis that land cover is not likely to be an important driver of yield mediated by either ecosystem services or disservices, and there are other factors that are far more influential predictors of yield. Among these, development of better crop varieties has long been the target of research (Reynolds et al., 2012), though the presence of any appreciable variance in our results implies that farmers may routinely be making suboptimal choices of which variety to plant for their district. This variance may also be attributable to intrinsic differences in the seed weights of varieties marketed for different applications, and warrants further analysis (e.g., malting, food and feed barley; Baik and Ullrich, 2008). Year-to-year variability explains additional variance, and is likely primarily driven by seasonal differences in weather given the synchrony observed among crops (Fig. A1). Another important factor is spatial variability, which stands in for unmodelled abiotic and biotic environmental conditions, and is expected to be considerable across a large geographic region (Fig. A2).

landscape complexity and yields in multiple crops across a large agricultural region (n = 10,069 yield records covering 0.25 million km2 seeded area) where environmental conditions are highly variable and span a broad range of land use intensities from districts with high proportions of crop and others dominated by non-crop land covers (e.g., Fig. 2, C). It is, to our knowledge, the largest area over which this relationship has been identified. Even though the practical contribution of non-crop land cover to yield was small given the limited variance it explained (ICC < 0.01; Table 1), the geographic area over which this relationship holds implies that these findings may have important consequences for the total amount of crop harvested. In particular, the findings may have practical relevance for the conservation of non-crop land covers across North American temperate grasslands. Agricultural and land use practices are consistent throughout the Canadian prairies, and much of the variation in land cover, moisture, and vegetation across the region is represented in this analysis. A geographical extrapolation of these results to the crop rotation footprint of the Canadian prairie region (more than 500,000 km2; Agriculture and Agri-Food Canada, 2015) could therefore have continental significance. The potential to increase landscape complexity through the redesign of agricultural landscapes has been the topic of much attention (Bommarco et al., 2013; Landis, 2017). Relationships between landscape characteristics and the arthropods providing ecosystem services to agriculture have received considerable focus, and landscape interventions to return biodiversity have also been explored (Kleijn et al., 2009; Landis et al., 2000; Tscharntke et al., 2005). However, much more knowledge about the relationship between landscape and crop yield is likely needed to drive implementation at a meaningful geographic scale. Our study helps by comparing crop yield potential across a broad gradient of landscape complexity, however farmers will need more guidance. While much has been studied and modelled about what to plant and restore (e.g., Kremen and M’Gonigle, 2015; Lonsdorf et al., 2009), connections to observed crop yields are often missing. Future studies should seek to provide clear instruction on how yield is related to the amount and spatial configuration of land cover and relate this to the specific ecosystem services that cover provides. In addition, not all patches of non-crop land cover are created equal. In our study we made a simplifying assumption to treat all non-crop land cover equally, a necessary step because land cover products with appropriate spatial and thematic resolution do not exist for this large extent. A binary land cover variable may also be advantageous, as it offers interpretability at the large extent of this study where classes of non-crop land covers are predominantly vegetation, and the types of vegetation differ geographically. However, exploring the relationship between the type of non-crop land cover and crop productivity will improve our understanding of the mechanisms at work. At smaller scales than those we studied, site conditions are likely to influence the amount and the direction of the effect on agricultural production (e.g., Duriaux Chavarría et al., 2018). For example, a particular patch of non-crop land cover may support a suite of both ecosystem services and perhaps disservices, with the net effect on yield reflecting a balance particular to the biotic and abiotic conditions supported by that land cover. In highly mechanized agroecosystems, precision agriculture (Mulla, 2013) may be helpful for understanding such differences among patch types, where data from precision yield monitors can be used to spatially-localize influences of specific non-crop land covers for crop productivity. While we found evidence for the prediction that non-crop land

4.2. Significance The main result is novel because it identifies an association between

10

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

cover in a district is positively associated with crop yields in that district, it would be ecological fallacy (sensu Robinson, 2009) to extend relationships identified at the aggregate scale to the field level. Thus, we do not imply that our results instruct farmers to adopt off-field practices that introduce more non-crop land covers into or near their fields in order to improve yields, nor that such steps are certain to be free of pest, weed, moisture or shade disservices to crops. Rather, we view these findings as identifying a pattern in a large-scale observational setting. Further, caution is warranted in extrapolating this model beyond the range of land cover and landscape complexities that were used to build it. We were unable to access field-level data for this analysis due to privacy restrictions. However, acquiring data for individual fields will be essential to demonstrate whether the influence of non-crop land cover persists at the spatial scale that matters to individual farmers, and is of sufficient magnitude to warrant policies intended to promote it. Such data exist, and in Canada are held by crop insurance providers, as well as government agencies. An alternative is the enrollment of farmers in a voluntary program of data sharing, though this may present challenges to obtaining coverage at meaningful geographic extents.

cover within annual cropland may need to be leveraged (Asbjornsen et al., 2014; Scherr and McNeely, 2008). A key finding of our study is that increasing non-crop land cover is unlikely, on average, to have negative consequences for crop yield in the Canadian prairie temperate grassland region, and at best may improve yields per unit area. This suggests that by increasing landscape complexity, farmers have the potential to receive a tangible yield credit for their land conservation contributions, while harnessing ecosystem services, creating habitat for biodiversity, and maintaining perennial carbon stores. Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements Agriculture Financial Services Corporation (AFSC) provided yield data. H. Bloom and M. Canuel assisted with land cover verification. Funding for this work was provided by the Alberta Canola Producers Commission, Manitoba Canola Growers Association, SaskCanola, Ducks Unlimited Canada’s Institute for Wetland and Waterfowl Research and Natural Sciences and Engineering Research Council (NSERC).

5. Conclusion In the context of global change, even small incentives for the establishment, retention or restoration of perennial and semi-natural land Appendix A .

Fig. A1. Inter-annual variation in yield over-plotted from separate regressions for each of seven crops. These marginal smooths suggest that among-season variation in yield may be synchronized across crops.

11

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Fig. A2. Marginal Markov random field smooths from separate regressions on each of seven crops. These smooths model covariation among adjacent municipal districts and when mapped likely demonstrate patterns in environmental variables such as soil and climate that are influencing yield.

12

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Table A1 Spatial data inputs used to produce a binary land cover raster map for Alberta, Canada, distinguishing crop and non-crop land covers at 30 m spatial resolution. In all cases spatial resolution was equal to, or coarser than, the original spatial data sources. Raster was compiled by overlaying these spatial data sources in the order given, after aggregating land cover classes of similar types. Spatial data source

Features represented

Reference

Agriculture and Agri-Food Canada Crop Inventory (2017)

Crop land covers (e.g., fields by crop type; pasture and hay lands)

Alberta Biodiversity Monitoring Institute Human Footprint Layer (2016) Alberta Merged Wetland Inventory

Non-crop land covers (e.g., roads, rail corridors, urban, natural and disturbed vegetated surfaces) Non-crop land covers (wetlands)

OpenStreetMap (2017)

Non-crop land covers (roads; rail corridors)

Annual Space-Based Crop Inventory for Canada (2017). Centre for Agroclimate, Geomatics and Earth Observation, Science and Technology Branch, Agriculture and Agri-Food Canada. https://open.canada.ca/data/en/dataset/ba2645d54458-414d-b196-6303ac06c1c9 Wall-to-Wall Human Footprint Inventory. (2016). Alberta Biodiversity Monitoring Institute. http://abmi.ca/home/data-analytics/da-top/da-productoverview/GIS-Land-Surface/HF-inventory.html Alberta Merged Wetland Inventory. (2017). Alberta Environment and Parks, Government of Alberta. https://maps.alberta.ca/genesis/rest/services/Alberta_ Merged_Wetland_Inventory/Latest/MapServer/ OpenStreetMap contributors. (2017). https://planet.osm.org/

Table A2 Crop yield data sample sizes by crop, year, Alberta municipal district and variety. Data were reported to a crop insurance program and reflect a precise accounting of yield per seeded acre. Year CANOLA

BARLEY

PEA

OAT

WHEAT (Hard Red Spring)

WHEAT (Canadian Prairie Spring)

WHEAT (Soft White Spring)

2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total 2012 2013 2014 2015 2016 2017 Total

nMunicipal

nCrop

districts

variety

65 64 64 65 63 64

72 58 55 58 60 62

60 59 59 59 59 58

38 33 30 33 32 28

53 53 53 59 59 60

17 16 18 20 17 12

43 43 48 53 55 44

6 8 7 10 11 8

59

35

60 59 60 59 59

33 35 36 30 30

35

8

41 42 40 40 36

9 8 11 9 7

14

2

21 16 18 20 15

2 2 2 3 3

13

nTotal

Total crop (106 acres)

713 749 701 636 623 679 4101 387 403 335 353 353 280 2111 108 97 105 122 137 150 719 70 79 78 88 88 61 464 360

4.08 3.97 4.22 3.91 3.78 4.17 24.12 2.05 2.09 1.74 1.82 1.84 1.35 10.9 0.57 0.61 0.75 0.91 1.22 1.08 5.14 0.11 0.13 0.14 0.16 0.16 0.11 0.8 3.27

335 326 340 357 364 2082 57

3.43 3.19 3.16 2.84 2.92 18.81 0.34

88 67 75 97 68 452 20

0.59 0.52 0.52 0.56 0.4 2.93 0.07

28 22 22 28 20 140

0.14 0.11 0.09 0.13 0.08 0.61

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Table A3 Classification error rates for the binary land cover product representing crop and non-crop cover classes, obtained by randomly selecting ten 2 km x 2 km regions positioned in each of four scenes chosen to represent the geographic extent of Alberta croplands. Proportional area of region misclassified as Scene 1 2 3 4

Nregions

Crop mean

10 10 10 10

0.009 0.044 0.010 0.019

(se) 0.002 0.006 0.002 0.004

Non-crop mean

(se)

0.011 0.025 0.048 0.030

0.004 0.005 0.012 0.006

Table A4 Modelling of non-crop land cover variation across municipal districts. (a) Significant municipal district smooth terms (P < 0.001) imply that the function estimating continuous variation in non-production land cover differed from zero, where higher edf values indicate greater non-linearity. (b) Significant parametric terms (P < 0.001) imply that the average amount of non-production land cover in a municipal district (i.e. averaged over all annulus distances) differed from the reference municipal district mean. (a) Smooth terms edfa Intercept Northing × Easting (tensor product smooth) (average field size) (spline smooth) Municipal district mean (spline smooth) Municipal district effects (spline smooth and parametric) ACADIA ATHABASCA BARRHEAD BEAVER BIG LAKES BIRCH HILLS BONNYVILLE BRAZEAU CAMROSE CARDSTON CLEAR HILLS CLEARWATER CYPRESS FAIRVIEW FLAGSTAFF FOOTHILLS FORTY MILE GRANDE PRAIRIE GREENVIEW KNEEHILL LAC LA BICHE LAC STE. ANNE COUNTY LACOMBE LAMONT LEDUC LESSER SLAVE RIVER LETHBRIDGE MACKENZIE MINBURN MOUNTAIN VIEW NEWELL NORTHERN LIGHTS NORTHERN SUNRISE PAINTEARTH PARKLAND COUNTY PEACE PINCHER CREEK PONOKA PROVOST RED DEER

1669.25*

(b) Parametric terms Estimate −1.09*

1776.91b,* 3.38* 6.53*,c 0.03 11.2* 6.63* 11.76* 7.6* 19.57* 0.08 2.04 10.8* 9.09* 7.12* 9.92* 11.6* 11.9* 16.94* 11.88* 11.26* 13.35* 12.6* 7.35* 15.99* 13.62* 11.45* 0.07 5.91* 17.44* 10.05* 11.8* 13.76* 4.78* 14.1* 5.19* 14.55* 5.83* 9.37* 14.22* 10.16* 8.97* 7.25*

c

0.35* 0.19* 0.12* 0.40* 0.06* 0.35* 0.46* 0.00 0.18* 0.24* 0.35* 0.33* 0.03* 0.05* 0.18* 0.06* 0.19* 0.35* −0.24* 0.44* 0.37* 0.07* 0.00 0.14* 0.44* −0.15* 0.19* 0.03* 0.07* 0.28* 0.25* 0.13* 0.28* 0.32* 0.15* 0.28* 0.20* 0.25* 0.06*

(continued on next page)

14

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al.

Table A4 (continued) (a) Smooth terms edfa ROCKY VIEW SADDLE HILLS SMOKY LAKE SMOKY RIVER SPECIAL AREA 02 SPECIAL AREA 03 SPECIAL AREA 04 SPIRIT RIVER ST. PAUL STARLAND STETTLER STRATHCONA COUNTY STURGEON TABER THORHILD TWO HILLS VERMILION RIVER VULCAN WAINWRIGHT WARNER WESTLOCK WETASKIWIN WHEATLAND WILLOW CREEK WOODLANDS YELLOWHEAD

9.74* 12.8* 9.17* 14.94* 9.98* 0.31 7.58* 8.14* 19.07* 9.58* 6.42* 10.77* 2.91 13.44* 10.29* 11.35* 12.56* 3.26* 7.44* 10.78* 3.89* 4.34* 13.44* 9.64* 8.4* 5.78*

(b) Parametric terms Estimate 0.09* 0.23* 0.25* −0.13* 0.40* 0.36* 0.36* −0.05* 0.32* 0.02* 0.28* 0.27* 0.03* 0.00 0.20* 0.13* 0.15* −0.09* 0.19* 0.05* 0.07* 0.24* −0.07* 0.13* 0.40* 0.43*

* P < 0.001. a Estimated degrees of freedom. b Note that in a function-on-scalar regression a single variable smooth is entered as a tensor product (i.e. two-dimensional smooth), resulting in a much higher estimated degrees of freedom. c Reference municipal district for parametric effects.

References

Ekström, G., Ekbom, B., 2011. Pest control in agro-ecosystems: an ecological approach. Crit. Rev. Plant Sci. 30, 74–94. Fahrig, L., Baudry, J., Brotons, L., Burel, F.G., Crist, T.O., Fuller, R.J., Sirami, C., Siriwardena, G.M., Martin, J.-L., 2011. Functional landscape heterogeneity and animal biodiversity in agricultural landscapes. Ecol. Lett. 14, 101–112. Fischer, J., Abson, D.J., Butsic, V., Chappell, M.J., Ekroos, J., Hanspach, J., Kuemmerle, T., Smith, H.G., von Wehrden, H., 2014. Land sparing versus land sharing: moving forward. Conserv. Lett. 7, 149–157. Foley, J.A., Defries, R., Asner, G.P., Barford, C., Bonan, G., Carpenter, S.R., Chapin, F.S., Coe, M.T., Daily, G.C., Gibbs, H.K., Helkowski, J.H., Holloway, T., Howard, E.A., Kucharik, C.J., Monfreda, C., Patz, J.A., Prentice, I.C., Ramankutty, N., Snyder, P.K., 2005. Global consequences of land use. Science 309, 570–574. Forman, R.T.T., 2009. Urban ecology: Science of cities, Urban Ecology: Science of Cities. Gallé, R., Happe, A., Baillod, A.B., Tscharntke, T., Batáry, P., 2019. Landscape configuration, organic management, and within‐field position drive functional diversity of spiders and carabids. J. Appl. Ecol. 56, 63–72. Garbach, K., Milder, J.C., DeClerck, F.A.J., Montenegro de Wit, M., Driscoll, L., GemmillHerren, B., 2017. Examining multi-functionality for crop yield and ecosystem services in five systems of agroecological intensification. Int. J. Agric. Sustain. 15, 11–28. Garibaldi, La., Carvalheiro, L.G., Vaissière, B.E., Gemmill-herren, B., Hipólito, J., Freitas, B.M., Ngo, H.T., Azzu, N., Sáez, A., Åström, J., An, J., Blochtein, B., 2016. Mutually beneficial pollinator diversity and crop yield outcomes in small and large farms. Science 351, 388–391. Garibaldi, L.A., Steffan-Dewenter, I., Kremen, C., Morales, J.M., Bommarco, R., Cunningham, S.A., Carvalheiro, L.G., Chacoff, N.P., Dudenhöffer, J.H., Greenleaf, S.S., Holzschuh, A., Isaacs, R., Krewenka, K., Mandelik, Y., Mayfield, M.M., Morandin, L.A., Potts, S.G., Ricketts, T.H., Szentgyörgyi, H., Viana, B.F., Westphal, C., Winfree, R., Klein, A.M., 2011. Stability of pollination services decreases with isolation from natural areas despite honey bee visits. Ecol. Lett. 14, 1062–1072. Gene, S.M., Hoekstra, P.F., Hannam, C., White, M., Truman, C., Hanson, M.L., Prosser, R.S., 2019. The role of vegetated buffers in agriculture and their regulation across Canada and the United States. J. Environ. Manage. 243, 12–21. Godfray, H.C.J., Beddington, J.R., Crute, I.R., Haddad, L., Lawrence, D., Muir, J.F., Pretty, J., Robinson, S., Thomas, S.M., Toulmin, C., 2010. Food security: the challenge of feeding 9 billion people. Science 327, 812–818. Hijmans, R.J., 2015. Raster: Geographic Data Analysis and Modelling. R Package Version 2.4-20. http://CRAN.R-project.org/package=raster. Huang, L., Scheipl, F., Goldsmith, J., Gellar, J., Harezlak, J., McLean, M.W., Swihart, B., Xiao, L., Crainiceanu, C., Reiss, P., 2016. Refund: regression with functional data. R package version 0.1-14. Hungate, B.A., Barbier, E.B., Ando, A.W., Marks, S.P., Reich, P.B., van Gestel, N., Tilman, D., Knops, J.M.H., Hooper, D.U., Butterfield, B.J., 2017. The economic value of grassland species for carbon storage. Sci. Adv. 3, e1601880.

Agriculture and Agri-Food Canada, 2015. AAFC Annual Crop Inventory. Government of Canada. Alignier, A., Bretagnolle, V., Petit, S., 2012. Spatial patterns of weeds along a gradient of landscape complexity. Basic Appl. Ecol. 13, 328–337. Asbjornsen, H., Hernandez-Santana, V., Liebman, M., Bayala, J., Chen, J., Helmers, M., Ong, C.K., Schulte, L.A., 2014. Targeting perennial vegetation in agricultural landscapes for enhancing ecosystem services. Renew. Agric. Food Syst. 29, 101–125. Baik, B.-K., Ullrich, S.E., 2008. Barley for food: characteristics, improvement, and renewed interest. J. Cereal Sci. 48, 233–242. Benton, T.G., Vickery, J.A., Wilson, J.D., 2003. Farmland biodiversity: is habitat heterogeneity the key? Trends Ecol. Evol. (Amst.) 18, 182–188. Bjørnstad, O.N., 2016. Ncf: Spatial Nonparametric Covariance Functions. R Package Version. https://doi.org/10.1088/0953-8984/18/18/S02. Boesing, A.L., Nichols, E., Metzger, J.P., 2017. Effects of landscape structure on avianmediated insect pest control services: a review. Landsc. Ecol. 32, 931–944. Bommarco, R., Kleijn, D., Potts, S.G., 2013. Ecological intensification: harnessing ecosystem services for food security. Trends Ecol. Evol. (Amst.) 28, 230–238. Bommarco, R., Vico, G., Hallin, S., 2018. Exploiting ecosystem services in agriculture for increased food security. Glob. Food Sec. 17, 57–63. Chaplin‐Kramer, R., O’Rourke, M.E., Blitzer, E.J., Kremen, C., 2011. A meta‐analysis of crop pest and natural enemy response to landscape complexity. Ecol. Lett. 14, 922–932. Dalu, T., Wasserman, R.J., Dalu, M.T.B., 2017. Agricultural intensification and drought frequency increases may have landscape‐level consequences for ephemeral ecosystems. Glob. Chang. Biol. 23, 983–985. DigitalGlobe, 2016. World Imagery. Doherty, K.E., Howerter, D.W., Devries, J.H., Walker, J., 2018. Prairie pothole region of North America. In: Finlayson, C., Milton, G.R., Prentice, R.C., Davidson, N.C. (Eds.), The Wetland Book. Springer, Dordrecht, Netherlands, pp. 1–10. Duarte, G.T., Santos, P.M., Cornelissen, T.G., Ribeiro, M.C., Paglia, A.P., 2018. The effects of landscape patterns on ecosystem services: meta-analyses of landscape services. Landsc. Ecol. 33, 1247–1257. Duriaux Chavarría, J.Y., Baudron, F., Sunderland, T., 2018. Retaining forests within agricultural landscapes as a pathway to sustainable intensification: evidence from Southern Ethiopia. Agric. Ecosyst. Environ. 263, 41–52. Ecological Stratification Working Group, 1995. A National Ecological Framework for Canada. A National Ecological Framework for Canada. Egli, L., Meyer, C., Scherber, C., Kreft, H., Tscharntke, T., 2018. Winners and losers of national and global efforts to reconcile agricultural intensification and biodiversity conservation. Glob. Chang. Biol. 24, 2212–2228.

15

Agriculture, Ecosystems and Environment 290 (2020) 106724

P. Galpern, et al. Karp, D.S., Frishkoff, L.O., Echeverri, A., Zook, J., Juárez, P., Chan, K.M.A., 2018. Agriculture erases climate‐driven β‐diversity in Neotropical bird communities. Glob. Chang. Biol. 24, 338–349. Kleijn, D., Kohler, F., Báldi, A., Batáry, P., Concepción, E.D., Clough, Y., Díaz, M., Gabriel, D., Holzschuh, A., Knop, E., 2009. On the relationship between farmland biodiversity and land-use intensity in Europe. Proc. R. Soc. Lond. B Biol. Sci. 276, 903–909. Kokoszka, P., Reimherr, M., 2017. Introduction to Functional Data Analysis. CRC Press. Kremen, C., M’Gonigle, L.K., 2015. Small-scale restoration in intensive agricultural landscapes supports more specialized and less mobile pollinator species. J. Appl. Ecol. 52, 602–610. Lamb, A., Green, R., Bateman, I., Broadmeadow, M., Bruce, T., Burney, J., Carey, P., Chadwick, D., Crane, E., Field, R., 2016. The potential for land sparing to offset greenhouse gas emissions from agriculture. Nat. Clim. Chang. 6, 488. Landis, D.A., 2017. Designing agricultural landscapes for biodiversity-based ecosystem services. Basic Appl. Ecol. 18, 1–12. Landis, D.A., Wratten, S.D., Gurr, G.M., 2000. Habitat management to conserve natural enemies of arthropod pests in agriculture. Annu. Rev. Entomol. 45, 175–201. Larmour, J., 2005. Laying Down the Lines: a History of Land Surveying in Alberta. Brindle and Glass. Laurance, W.F., Sayer, J., Cassman, K.G., 2014. Agricultural expansion and its impacts on tropical nature. Trends Ecol. Evol. (Amst.) 29, 107–116. Letourneau, D.K., Armbrecht, I., Rivera, B.S., Lerma, J.M., Carmona, E.J., Daza, M.C., Escobar, S., Galindo, V., Gutiérrez, C., López, S.D., 2011. Does plant diversity benefit agroecosystems? A synthetic review. Ecol. Appl. 21, 9–21. Lonsdorf, E., Kremen, C., Ricketts, T., Winfree, R., Williams, N., Greenleaf, S., 2009. Modelling pollination services across agricultural landscapes. Ann. Bot. 103, 1589–1600. https://doi.org/10.1093/aob/mcp069. Marx, B.D., Eilers, P.H.C., 1999. Generalized linear regression on sampled signals and curves: a P-spline approach. Technometrics 41, 1–13. McGarigal, K., Cushman, S.A., Neel, M.C., Ene, E., 2012. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps. University of Massachusettes, Amherst, MA https://doi.org/citeulike-article-id:287784. Mulla, D.J., 2013. Twenty five years of remote sensing in precision agriculture: key advances and remaining knowledge gaps. Biosyst. Eng. 114, 358–371. O’Reilly-Nugent, A., Palit, R., Lopez-Aldana, A., Medina-Romero, M., Wandrag, E., Duncan, R.P., 2016. Landscape effects on the spread of invasive species. Curr. Landsc. Ecol. Rep. 1, 107–114. Phalan, B., Onial, M., Balmford, A., Green, R.E., 2011. Reconciling food production and biodiversity conservation: land sharing and land sparing compared. Science 333, 1289–1291. Ponisio, L.C., M’gonigle, L.K., Kremen, C., 2016. On‐farm habitat restoration counters biotic homogenization in intensively managed agriculture. Glob. Chang. Biol. 22, 704–715. Reynolds, M., Foulkes, J., Furbank, R., Griffiths, S., King, J., Murchie, E., Parry, M., Slafer, G., 2012. Achieving yield gains in wheat. Plant Cell Environ. 35, 1799–1823.

Robinson, W.S., 2009. Ecological correlations and the behavior of individuals. Int. J. Epidemiol. 38, 337–341. Rusch, A., Chaplin-Kramer, R., Gardiner, M.M., Hawro, V., Holland, J., Landis, D., Thies, C., Tscharntke, T., Weisser, W.W., Winqvist, C., 2016. Agricultural landscape simplification reduces natural pest control: a quantitative synthesis. Agric. Ecosyst. Environ. 221, 198–204. Scherr, S.J., McNeely, J.A., 2008. Biodiversity conservation and agricultural sustainability: towards a new paradigm of “ecoagriculture” landscapes. Philos. Trans. Biol. Sci. Smith, P., 2014. Do grasslands act as a perpetual sink for carbon? Glob. Chang. Biol. 20, 2708–2711. Statistics Canada, 2018. Table 32-10-0359-01 Estimated Areas, Yield, Production, Average Farm Price and Total Farm Value of Principal Field Crops, in Metric and Imperial Units [WWW Document]. https://www150.statcan.gc.ca/t1/tbl1/en/tv. action?pid=3210035901. Titeux, N., Henle, K., Mihoub, J., Regos, A., Geijzendorffer, I.R., Cramer, W., Verburg, P.H., Brotons, L., 2016. Biodiversity scenarios neglect future land‐use changes. Glob. Chang. Biol. 22, 2505–2515. Tscharntke, T., Clough, Y., Wanger, T.C., Jackson, L., Motzke, I., Perfecto, I., Vandermeer, J., Whitbread, A., 2012. Global food security, biodiversity conservation and the future of agricultural intensification. Biol. Conserv. 151, 53–59. Tscharntke, T., Klein, A.M., Kruess, A., Steffan‐Dewenter, I., Thies, C., 2005. Landscape perspectives on agricultural intensification and biodiversity–ecosystem service management. Ecol. Lett. 8, 857–874. Tubiello, F.N., Salvatore, M., Ferrara, A.F., House, J., Federici, S., Rossi, S., Biancalani, R., Condor Golec, R.D., Jacobs, H., Flammini, A., 2015. The contribution of agriculture, forestry and other land use activities to global warming, 1990–2012. Glob. Chang. Biol. 21, 2655–2660. Vickruck, J.L., Best, L.R., Gavin, M.P., Devries, J.H., Galpern, P., 2019. Pothole wetlands provide reservoir habitat for native bees in prairie croplands. Biol. Conserv. 232, 42–50. White, E.V., Roy, D.P., 2015. A contemporary decennial examination of changing agricultural field sizes using Landsat time series data. Geo Geogr. Environ. https://doi. org/10.1002/geo2.4. Wickham, H., 2016. ggplot2 – Elegant Graphics for Data Analysis (use R), Elegant Graphics for Data Analysis. Srpinger-Verlag, New York. Williams, D.R., Phalan, B., Feniuk, C., Green, R.E., Balmford, A., 2018. Carbon storage and land-use strategies in agricultural landscapes across three continents. Curr. Biol. 28, 2500–2505. Winqvist, C., Bengtsson, J., Aavik, T., Berendse, F., Clement, L.W., Eggers, S., Fischer, C., Flohre, A., Geiger, F., Liira, J., 2011. Mixed effects of organic farming and landscape complexity on farmland biodiversity and biological control potential across Europe. J. Appl. Ecol. 48, 570–579. Wood, S.N., 2017. Generalized Additive Models: an Introduction With R. Chapman and Hall/CRC.

16