Agro-biodiversity has increased over a 95 year period at sub-regional and regional scales in southern Quebec, Canada

Decline in agricultural biodiversity (cultivated species and wild species used for food or that support agro-ecosystem functioning) at the farm scale has fueled concerns about potential negative effects of this biodiversity loss on the ecological and economic sustainability of agro-ecosystems. Despite these concerns, formal assessment of how agro-biodiversity has historically changed at scales larger than individual farms is fragmented. We quantified the changes in the abundance of 10 crop and livestock species, their overall diversity, and the way they were mixed in ‘baskets’ of agricultural products from 1911 to 2006, at a sub-regional (15 Regional County Municipalities) and regional scales. We found that the diversity of agricultural products increased at the regional scale. From 1911 to 1966, the region produced fodder, milk and maple, mixed in two low-diversity baskets. After 1966, the region provided a greater variety of baskets composed of newly introduced cash crops and high-value livestock. All baskets provided were themselves more diverse than historically and varied greatly in composition across space. Increasing regional diversity was related to changes in agricultural policy, while the variation in the composition of baskets produced was related to biophysical and socioeconomic characteristics. Our results indicate that agricultural transformations of the 21th century did not invariably lead to agro-biodiversity loss at large scales. We have also demonstrated that combining diversity measures at multiple scales with the analysis of compositional change of agricultural products over long time periods could improve research on the links between agro-biodiversity dynamics and resilience.


Introduction
Agro-biodiversity, which includes cultivated species and wild species used for food or that support agroecosystem functioning (Frison et al 2011), has declined at the farm scale in recent decades, primarily due to increasing agricultural industrialization (Matson et al 1997, Tilman 2002, Hammer and Khoshbakht 2005, Tscharntke et al 2005, Frison et al 2011. With advances in mechanization and increasing availability of chemical inputs over the last 50 years, it became more profitable for farmers to produce a limited set of crop or livestock species, often managed in monocultures with low genetic diversity (Matson et al 1997, Swift et al 2004. The divorce of agriculture from the ecological processes driven by agro-biodiversity (including nutrient and water cycling, pest and disease control, pollination) has fueled concerns about the long-term sustainability of production systems. These concerns relate to the capacity of agro-ecosystems to sustain high and stable yields (Tilman 2002), support a diversified diet (Frison et al 2011, Khoury et al 2014, cope with economic shocks and sustain livelihoods (Rotz and Fraser 2015).
Much empirical evidences supports on-farm 'diversification' as a sustainability principle to counter the potential negative effects of agro-biodiversity loss (see Lin 2011 for a review). Manipulative experiments have showed that fields with multiple genotypes of one crop species could increase yield, suppress disease (Zhu et al 2000), and stabilize producers' income (Di Falco and Chavas 2006). Systems with higher crop species diversity also support greater wild biodiversity, generate multiple ecosystem services (Kremen andMiles 2012, McDaniel et al 2014), are more productive (Smith et al 2008, Davis et al 2012 and more resilient to climate change than monocultures (Snapp et al 2010, Kremen and Miles 2012, Gaudin et al 2015. Despite interests in the effects of on-farm agrobiodiversity loss and research on ways to counter them, few studies have quantitatively assessed changes in agro-biodiversity across multiple spatial scales (Swift et al 2004). Among these few, two recent studies revealed that the direction of historical trends depend on the spatial scale. Khoury et al (2014) showed that countries globally increased the diversity of crop species contributing to national food supplies between 1961 and 2009. In the USA, Aguilar et al (2015) found a decline in crop species diversity at the national scale between 1978 and 2012, but revealed both declining and increasing diversity trends among production regions, and among counties within regions. Further examination of how historical transformations of agriculture may have affected agro-biodiversity at different spatial scales is critical, especially to understand the implications of past trends for agricultural resilience. Indeed, as scale widens, natural and anthropogenic disturbances become more frequent and intense, and more difficult and costly to manage (Swift et al 2004). The role of agro-biodiversity in supporting resilience and adaptation of agricultural systems in the face of these stressors might thus also heighten with scale.
We quantified changes in the production of ten crop and livestock species over a 95 year period at county and regional scales in southern Quebec, Canada. Our work specifically addresses the following questions: (1) How did agro-biodiversity-including individual species abundance, and the mix of crop and livestock produced together-change through time from 1911-2006 in this region, (2) How did the trajectories of change in agro-biodiversity vary across space in this region during this time period? and (3) How did local biophysical and socioeconomic factors relate to changes in agricultural production? Assessing these factors is critical to understanding the spatio-temporal dynamics of agro-ecosystems, which are fundamentally one of the most closely coupled human and natural systems (Tomich et al 2011).

Material and methods
In this study, we use the term 'agricultural products' to refer to the ten crop, livestock species and maple products quantified. The term 'basket' refers to the mix of agricultural products produced together at the same place and at the same time. For example, a basket characteristic of a dairy region could comprise different fodder crops like hay and oats, as well as milk cows. Finally, we assessed how the abundance of agricultural products and their diversity within baskets changed using government census records across 15 administrative units in a region of southern Quebec from 1911 to 2006.

Study area
Our study area is the Monteregie (11 111 km 2 ), a region of southern Quebec, Canada, located east of Montreal. The region is divided into 15 administrative units, called Regional County Municipalities (RCMs, figure 1). It belongs to the humid-cold temperate ecoclimatic region of Canada (EcoRegions Working Group 1989), and its vegetation is part of the northern hardwoods forest type (Boucher et al 2006). The topography of the region is relatively flat (ranges from ∼200 to 400 m above sea level), except for the series of seven Monteregian hills, formed of intrusive igneous rocks, which are somewhat higher. The lowland plains consist of nutrient-rich clay deposits formed by the post-glacial Champlain Sea.
The Monteregie has undergone several major land transformations since European settlement in the early 1800s. After an intensive period of logging activity, starting around 1870, the landscape was converted into an agricultural matrix with only small remnant woodlots Bouchard 2007, Simard andBouchard 1996). The Monteregie of today is characterized by a mixed landscape dominated by agricultural lands (75%), urban centers (10%) and remnant forest areas (15%) (Létourneau et al 2009).
References of data sets and maps used for reconstructing agricultural changes and assess socioeconomic and biophysical variables are detailed in the supplementary online material.
Reconstructing long-term agricultural changes using census records We used the Canadian Agricultural Census to reconstruct agricultural production within RCMs (referred to as census divisions in the census) across the Monteregie. We collected data on the production of five crops (oats, hay, maize, soybeans, and wheat), three livestock species used for meat (pork, horned cattle and poultry which comprises hens, chicken, duck, turkey and geese), milk cows (used to estimate milk production), and maple syrup (table 1). Based on the census data, we estimated that the crops we studied represent 64%-84% of the total cultivated area in the region, depending on the year. For animal production, the importance of dairy and pork production in the Monteregie has been highlighted in the literature (Ruiz and Domon 2009). Agricultural data from 1911 to 1951 were only available at 10 year intervals in hardcopy format and had to be entered manually. From 1961 to 2006, the data were available digitally at 5 year intervals.
Because RCM administrative boundaries have changed through time, we standardized our dataset using the 2006 digital boundary file as a reference. To do this, the boundary maps from 1911 to 2001 were individually overlaid with the 2006 reference map, and the overlapping boundaries were clipped to create polygons representing the smallest units of discrete area. We then re-aggregated the census divisions by reassigning each polygon and the proportional value of the agricultural products associated with its area, to match the 2006 boundaries. Because no digital boundary maps were available from 1961 to 1981, we created our own files by comparing the map of 2006 to the digital map of 1951 and editing the 1986 map according to the boundary changes for the missing census years.

Explaining agricultural changes by socioeconomic and biophysical variables
We used data on population density, distance to the main economic center (Montreal), and the capacity of the soil to support agriculture as explanatory variables to assess the spatiotemporal changes in agricultural production. The Canadian Census of Population provided population data at the same spatial and temporal resolution as for agriculture. We used a national inventory to estimate the soil potential and limitations for agriculture. This inventory defines eight soil classes (from 0 to 7): class 0 refers to organic soils, classes 1 to 4 to soils capable of sustained use for field crops with increasing limitations from class 1 to 4. Classes 5-7 are soils with severe limitations, or not capable of crop production. For each RCM, we determined the class of soil that covered the most area and assigned that value to the RCM. Finally, we calculated the linear distance (in km) from the geometric center of each RCM to Montreal using digital boundary shapefiles for 2006.

Statistical analysis
We replaced missing data (1.7% of the dataset) at the beginning (1911) or end (2006) of the time series with the most recent value for the same RCM, i.e., from 1921 and 2001 respectively, and used linear interpolation for all other dates. We used the x′=sqrt(x) transformation to meet assumptions of normality, and  standardized the data set to unit variance and zeromean 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 (R Development Core Team 2011). First, we examined historical changes in agricultural production at the regional scale. We used a principal component analysis (PCA) to visualize the production of all 10 products across all RCMs over the entire time period. We then examined how the regional item diversity changed over time by using the Shannon diversity index (H′). This index takes into account both the richness and abundance of products without favouring rare or common species. The Shannon index was calculated for each year based on the sum of each of the 10 agricultural products per year, and compared across RCMs. The Shannon index values were subsequently transformed into an effective number of agricultural products (ENAPs), computed following Jost (2006): ENAP=exp H′ . The resulting value is an estimate of the number of products dominating agricultural production in a particular RCM at a particular time step. The resulting diversity index facilitates cross-study comparisons as the values have the same unit and the same properties, no matter the diversity index it has been calculated from (Shannon index here).
Second, we examined historical changes in agricultural production at the RCM scale using a timeconstrained multivariate regression tree (MRT) (Legendre and Gauthier 2014). This analysis allowed us to identify (i) RCMs following the same trajectories of change in the basket of agricultural products produced, and (ii) the socioeconomic and biophysical factors explaining these patterns. The analysis splits the RCMs into clusters based on a matrix of explanatory variables. We used time (census year), distance from Montreal, population density, and soil agricultural capability as explanatory variables. The MRT results are represented graphically as a tree with discriminating explanatory variables identified at each node. Each final cluster represents a group of RCMs that produce a similar mix of agricultural products under similar conditions (i.e. time period, population density, soil quality and distance to Montreal). When a node contributed to less than 2% of the total variance, clusters resulting from that node were aggregated. Finally, we calculated the diversity of the mix of agricultural products provided by RCM in each cluster, per year, using the procedure detailed above (H′ transformed into ENAP).
While the MRT provides insights into similarities between RCMs in terms of the mix of agricultural products produced, we further examined how the diversity of products changed through time among individual RCMs (i.e. beta diversity). According to (Legendre et al 2010), beta diversity can be conceived as the total variation in the agricultural data and can be decomposed into the contribution of each RCM (the 'local contribution to beta diversity' (LCBD)). The LCBD indexes obtained are comparative indicators of the ecological uniqueness of each RCM in terms of the composition of the baskets. LCBD values were calculated using the Euclidean method and tested for significance by random, independent permutations of the columns of the agricultural products matrix. The null hypothesis (H 0 ) that is tested posits that agricultural products are distributed randomly, independently of one another among sites. An RCM characterized by a significantly high LCBD value provides a basket of agricultural products that is unique compared to the other sites. We conducted a quadratic linear model on log-transformed LCBD values to test whether these values changed significantly through time.

Results
The PCA revealed a shift in agricultural products over time at the regional scale (figures 2(a) and (b), 71% of variance captured on the two first axis). The production of hay, oats and milk, dominated all RCMs from 1911 to 1961, and were then replaced by wheat, soybeans and maize from 1966 to 2006. In addition to the shift in crops (shown along axis 1 in figure 2(b)), the production of livestock (shown on axis 2) increased across the region starting in 1966, fuelling an increasing separation between RCMs producing large amounts of livestock and those which produce less or none. At the regional scale, agricultural production diversified through time until 1991. The diversity of agricultural products across the study region was the lowest, and varied little among RCMs, from 1911 to 1961, with an average (±SD) values of 6.35±0.14 (figure 2(c)). The major increase in diversity occurred between 1966 and 1991, when the diversity reached its maximum of 8.79. Then, diversity decreased slightly from 1991 to 2006 (ENAP 2006 =8).
When looking at the combination of agricultural products produced by each MRC, we found that the diversity of baskets among RCMs (i.e., beta-diversity) significantly increased through time ( During the 65 year period from 1911 to 1976, the RCMs were grouped into only three clusters (clusters 1-3 ( figure 4(a))), indicating that most RCMs produced very similar baskets of agricultural products. After 1976, the RCMs provided more differentiated baskets of products and were grouped into six clusters (clusters 4-9, figure 4(a)).
The diversity of agricultural products within baskets also increased through time. The clusters of RCMs that represented the production landscape from 1911 to 1966 (clusters 1 and 2) provided among the least diversified baskets of agricultural items of all times (figures 4(b) and (d)). Both baskets were composed primarily of oats and hay, with cluster 2 additionally having a higher production of milk and maple syrup. This cluster comprised all the RCMs except the closest ones to Montreal (N=2, distance<29 km, figure 4(c)).
Between 1971 and 1976, all RCMs belonged to cluster 3, which was characterized by the most diversified basket of agricultural products (figures 4(b) and (d)). This basket was composed of both products characteristic of the 1911-1966 period (oats, hay, milk and maple) and of those becoming dominant in the 1981-2005 period (including wheat, cattle and poultry, figure 4(b)).
Diversity within baskets changed following different trajectories after 1981. After 1981, the distance to Montreal explained most of the differentiation among clusters of RCMs. RCMs closer to Montreal (N=7, distance<48 km, figure 4(c)) comprised cluster 4 from 1981-1991, and cluster 5 from 1996 to 2006 ( figure 4(c)). Although both clusters specialized in crop production, the identity of the main crop shifted from wheat in cluster 4 to soybeans in cluster 5. This shift was accompanied by a decrease in the diversity of the baskets, from, on average 7.23±0.63 (cluster 4) to 5.91±0.68 (cluster 5). The diversity of the last basket continued to decline from 1996 to 2006 (figure 1SI). The RCMs located further away from Montreal  (N=8, distance>48 km, figure 4(b)) were grouped in clusters 6-9 ( figure 4(a)). RCMs where soils are the least suited for agriculture also produced highly diversified baskets of agricultural products (cluster 6, ENAP=8.13±0.28; cluster 7, ENAP=7.80± 0.14, figure 4(d)). RCMs in cluster 6, between 48 and 79 km from Montreal, produced a mix of crops and livestock. RCMs located the furthest away from Montreal (distance>79 km) belonged to cluster 7, and produced a mix of maple syrup and livestock. Where soils were the best suited for agriculture (clusters 8 and 9), population density differentiated RCMs between those producing mainly crops (cluster 8, population density>71 hab km −2 ), and those producing both livestock and crops (cluster 9, population density<71 hab km −2 ). Despite differences in the composition of the baskets provided, clusters 7, 8 and 9 showed relatively similar and high diversity values (figures 4(b) and (d)).

Discussion
Our work provides empirical evidence of agricultural diversification over time and at multiple spatial scales in the region around Montreal, QC, Canada. Historically, low-diversity baskets of agricultural products were provided in RCMs across the entire study region. After 1976, the region provided a greater variety of baskets, and the baskets provided were themselves diverse. The spatial arrangement of baskets on the landscape was related to socio-economic and biophysical characteristics of the region. These results indicate that agricultural transformations were not always accompanied by loss of agro-biodiversity, and show the value of formally assessing spatiotemporal changes beyond the farm scale and over long time periods.
Our results reveal a story of increasing diversity in agricultural products at the regional scale. Specialization of agricultural production in Quebec from 1911 Threshold values are given for both branches at each node. The thickness of each branch is proportional to its contribution to the explained variance. (b) Baskets of agricultural products, as represented by flower diagrams, provided by each of the nine clusters of RCMs. Each of the ten products is a color-specific petal which is proportional to the relative mean abundance of the product within each basket (normalized to range between 0 and 1). to 1966 was supported by the federal government, which subsidized food transport of grains from the western to the eastern provinces, where livestock farming was concentrated (Châtillon 1976, Morisset 1987. In the 1960s, the Quebec government promoted the diversification of agriculture to enhance the competitiveness and productivity of the province and to achieve self-sufficiency (Domon et al 1993, Brouillette-Paradis 2010. We demonstrated that the diversification was underlined by a shift in the type of agricultural products, from hay and oats to cash crops and high-value livestock. Due to mechanization, the shift in the type of products has been accompanied by an increase in individual field size and building area in some parts of the region (Pan et al 1999, Ruiz and Domon 2009, Jobin et al 2014. Our study contributes to showing that agricultural modernization over the 21th century could have led to a variety of trajectories that cannot be assumed to be towards decline without formal quantification.
Diversification at the regional scale did not occur in the same way in every RCM. From a single basket of agricultural products produced across the region from 1971 to 1981, certain RCMs differentiated and produced a unique basket of products compared to other RCMs. By reducing the overall degree of product overlap between RCMs, such differentiation could have conferred an economic, competitive, advantage (Bradshaw 2004). Ecological theory also posits that larger compositional dissimilarity among local species communities could stabilize regional ecosystem properties through spatial asynchrony effects (Loreau et al 2003, Wang and Loreau 2014, 2016. Our results thus suggest that the trajectory of diversification and differentiation among RCMs may have conferred increased resilience across the region. Although the overall beta-diversity increased among RCMs through time, the diversity of three baskets (produced by clusters 4, 5, 8) decreased and became more similar in composition after 1981. Because RCMs producing these clusters are spatially close, their convergent trajectories toward a soybeandominated production could facilitate the spread of pests and pathogens (Margosian et al 2009), thus affecting the use of insecticides (Meehan et al 2011) at larger spatial scales, and increasing regional vulnerability to price volatility (Di Falco andChavas 2006, Abson et al 2013). But soybeans specialization affected RCMs with the highest quality soils (Domon and Bouchard 2007) and did not expand throughout the entire region. RCMs located further away from Montreal, with the least suitable soil conditions, supported the highest mixtures of crop or livestock, and maple syrup stands. These results indicate that heterogeneity in biophysical conditions and gradients in market access could be important factors for the maintenance of agro-biodiversity at the regional scale.
In rural areas, where people chose to settle might also be an overlaying driver, in addition to soil quality, in determining the composition of baskets. In RCMs with similarly good soil conditions, crop production was spatially associated with the highest population density, and livestock production with the lowest population density. This pattern could result from a spatial incompatibility between human communities and industrial animal production (suggested in Raudsepp-Hearne et al 2010, Renard et al 2015). Considering that rural areas are becoming bedroom communities for a growing number of urban workers in Quebec (Paquette and Domon 2001), the relationship between human communities and agriculture might become key in the future of agricultural landscapes.
While it is beyond the scope of this study to provide guidelines for tracking changes in the diversity of agricultural products, our approach raised two methodological points that should be considered in future research. The first point concerns the measure of diversity itself. The Shannon index has been a popular diversity index in the ecological literature, including in studies examining agro-biodiversity (Finke and Swanson 1973, Brown and Schulte 2011, Aguilar et al 2015. The Shannon index takes into account both the richness and relative abundance of each product, but is insensitive to composition, e.g. the specific identity of the products grown or raised. In our study, complementing the regional assessment of diversity with a PCA revealed that agricultural diversification was underpinned by the substitution of crop species and livestock types over time. Thus, insights on mechanisms of agricultural transformation may only be fully gained by combining diversity measures with the analysis of the dynamics of agricultural products.
Our study also demonstrated the limit of using a diversity index alone to compare agricultural production across space. Although RCMs in clusters 7 and 9 provided similarly diverse baskets of products (values of 7.80 and 7.81, respectively), their composition varied greatly. Mixing livestock and maple syrup production (cluster 7) can promote synergies between carbon storage, recreational activities (sugar shack and hunting), and food production (Renard et al 2015). Forest stands can also provide habitat for biodiversity, including natural enemies of crop pests (see a review by Power 2010). The crop/livestock mixture (i.e., cluster 9) on the other hand, could achieve high productivity and optimal nutrient recycling (Pearson and Ison 1987). How agricultural products are combined and distributed is critical with respect to the purported environmental benefits of diversity in agriculture (Altieri 1998, Gliessman 2007, Margosian et al 2009. Representing the baskets of products and assessing how their composition changes provides more complete information to relate agro-biodiversity change with resilience.
Finally, the time scale at which changes in agricultural systems are examined matters. Although time-consuming, acquiring census data spanning almost a century allowed us to fully capture the period of main changes, provide a past context to account for the magnitude, direction and determinants of such changes, and detect potential warning trends in the recent period. The most recent literature on history of agriculture rarely explores trends prior to the 1960s (Bradshaw 2004, Daniel 2005, Brown and Schulte 2011, Abson et al 2013, Aguilar et al 2015, probably limited by the availability of older data in digital format. Some of this work does not even extend more than a decade or two (Gillmor 1987, Lyons 1988, Edwards 1991. Older literature focused mainly on the 1940s-1970s period, missing the most recent developments in agricultural production (Finke and Swanson 1973, Winsberg 1980, 1982.

Conclusions
Our study contributes to a growing line of research quantifying the spatiotemporal diversity of agricultural products (Khoury et al 2014, Aguilar et al 2015. Studies guiding agricultural policy have mainly focused at a single spatial scale (typically the field or the farm, Swift et al 2004, Gabriel et al 2006, or spanned limited time frames. We showed the importance of providing both larger spatial and longer temporal context to fully capture agro-biodiversity changes. Indeed, our results revealed that agricultural products diversified regionally over time, and that the spatial differentiation at the sub-regional level that increased through time could provide a measure of resilience. We have also demonstrated that combining diversity measures at multiple scales with an analysis of compositional change of agricultural products could improve research on the links between agrobiodiversity dynamics and resilience.