Evaluating the Relationships of Inter-Annual Farmland Vegetation Dynamics with Biodiversity Using Multi-Spatial and Multi-Temporal Remote Sensing Data

: Agricultural landscapes are highly dynamic ecosystems, but the e ﬀ ects of temporal farmland vegetation dynamics on species diversity have not been widely studied. In 93 sample farm landscapes in eastern Ontario, Canada, biodiversity data for seven taxa were collected in 2011 and 2012, prior to the initiation of this study. The goal of this study was to determine if trends and variability in vegetation productivity detected in these sample landscapes using long-term archived moderate and coarse resolution remote sensing time series data are related to the measured biodiversity. Mid-summer Moderate Resolution Imaging Spectroradiometer (MODIS) (2000–2011) and Landsat 5 (1985–2011) Normalized Di ﬀ erence Vegetation Index (NDVI) data were used with the Thiel–Sen slope and Contextual Mann–Kendall trend analysis to identify pixels showing signiﬁcant trends. NDVI temporal metrics included 1) the percentage of pixels in each landscape with a signiﬁcant negative or positive trend, and 2) the temporal coe ﬃ cient of variation (CV) of both the mean and spatial CV of landscape NDVI. Larger areas of signiﬁcant positive NDVI trends were found in the sample landscapes than negative trends, the former being associated with agricultural intensiﬁcation or crop changes and the latter with smaller areas of natural vegetation removal. Landsat better-detected changes in individual ﬁelds or small areas of natural vegetation due to its much smaller pixel size. In addition, the longer Landsat time series showed a change in the NDVI trend from positive (1985–2000) to negative or a leveling o ﬀ (2000–2011) for many pixels. In biodiversity modeling, the Landsat temporal CV of NDVI was negatively correlated with 2011–2012 plant and beetle diversity, while plant biodiversity was positively correlated with the percentage of pixels in a sample landscape showing a signiﬁcantly positive NDVI trend. No signiﬁcant relationships were found using the MODIS data. This study shows that temporal trends and variability in farmland vegetation density derived from Landsat data are related to biodiversity for certain taxa and that such relationships should be considered along with the more commonly studied spatial landscape attributes in evaluating landscape-level impacts of farming on biodiversity.


Introduction
Biodiversity has declined significantly in recent decades, mainly due to habitat loss and habitat fragmentation [1]. Land use change, especially the transformation of natural and semi-natural lands to farmlands to provide food for increasing human populations, is one of the major drivers of global biodiversity loss. In Canada, the area of natural lands decreased from over 300 million ha in the 1950s to 240 million ha in the 2000s, while the area of agricultural lands increased from 24.3 million hectares in the 1950s to 37.8 million hectares in 2016 [2]. In addition, the intensification of farming in existing agricultural lands presents increased risks to biodiversity and species at risk and negatively affects natural habitats and ecosystems [3]. Agriculture intensification involves processes such as reduction or removal of natural vegetation at field margins, use of similar crop types in adjacent fields, increased frequency of intensive cropping or harvesting, and increased use of pest and weed control, fertilization, irrigation, and artificial drainage systems [4]. The capacity of Canadian farmlands in carrying wildlife species was relatively stable from 1986 to 1996; however, it declined from 1996 to 2011, with agricultural intensification being a leading cause of this decrease [5]. Recently, conservation and alternative agricultural practices are being adopted in some areas to revitalize the landscape, enhance biodiversity, and minimize damage to living organisms. At the farm level, organic farming, crop rotation, no-tillage systems, and management of natural field margins are some practices that promote biodiversity [6][7][8]. At the landscape level, which is the focal scale of this study, increasing spatial heterogeneity is often suggested as a means to enhance species richness and abundance [9,10].
Landscape heterogeneity in agricultural environments can be characterized by the spatial and temporal variation in land cover. Increased spatial heterogeneity in farmlands is generally characterized by smaller fields, diverse crop types, and the presence of natural and semi-natural habitat patches such as forests and hedgerows. Previous studies have demonstrated the effects of spatial heterogeneity in agricultural lands on the richness and abundance of farmland species [9,[11][12][13][14][15][16][17][18][19][20]. Temporal heterogeneity in agricultural landscapes can be represented by the phenological cycles of crops and other vegetation, as well as inter-annual vegetation dynamics due to crop changes, farm management practices, or weather variations. Several studies have evaluated the impacts of temporal landscape dynamics on biodiversity using either past land cover data and simulations of future landscapes [21][22][23].
The application of remote sensing in biodiversity modfling has been reviewed in several papers and books (e.g., [24][25][26][27][28][29]). The most common approach to derive landscape data for biodiversity studies has been to use land cover maps generated from automated or manual classification of remote sensing imagery, combined with field observations of land cover. For example, in a study related to this paper (same landscapes and biodiversity data), Fahrig et al. [20] evaluated the relationships of alpha, beta, and gamma diversity of farmland species (plants, birds, beetles, spiders, bees, butterflies, and syrphids) with farmland compositional heterogeneity (the Shannon crop type diversity index) and configurational heterogeneity (field size) derived from air photos and field observations in farmlands of eastern Ontario. Mean crop field size had a stronger effect on biodiversity than the Shannon diversity index, and the relationship was consistently negative for all the taxa, i.e., increased configurational heterogeneity (smaller field size) was associated with higher biodiversity. Wilson et al. [30], working in southern Ontario and Quebec, using Agriculture and Agri-Food Canada (AAFC) agricultural inventory land cover maps and other data, found that specialist forest and wetland bird species richness declined with increasing area of arable crops while grassland bird species richness increased with increasing area of both forage and arable crops. Thematic classification of land use and land cover (LULC) is, however, subject to error, which can impact metrics derived from classified maps. In temporal analysis of LULC change by comparing classified maps for two dates, significant errors may result. Extending this to the classification of several to many dates of time series data can result in a significant error propagation in LULC change or trend analysis.
An alternative approach to generating landscape data for biodiversity studies is to use remotely sensed reflectance data, including combinations of spectral band reflectance such as vegetation and moisture indices. The use of reflectance data is advantageous over the use of LULC maps as it avoids the additional classification step and the associated potential for error propagation. However, reflectance data are limited in that they can only represent generalized classes such as vegetation through the use of specific spectral bands or vegetation indices. Despite this limitation, reflectance data are often well-calibrated and more suitable to time series analysis, as in this study, than the use of a series of LULC maps. The spectral variation hypothesis [28,31] suggests that spatial heterogeneity (i.e., pixel-to-pixel variability) of the reflectance in the spectral bands or derived indices that are indicative of vegetation is associated with the spatial heterogeneity of vegetation, which is often related to habitat heterogeneity and biodiversity. For example, Normalized Difference Vegetation Index (NDVI) [(NIR − R) / (NIR + R)] is often used as a metric representing above ground primary productivity over a landscape and its spatial variability has been associated with biodiversity and species distribution [21,32]. In the research related to this study (the same biodiversity data but for one-half the landscapes of this study), Duro et al. [33] used Landsat NDVI and three NDVI local spatial variance metrics (Coefficient of Variation, Geary's C, and Moran's I) in modeling the diversity of birds, butterflies, and plants in eastern Ontario farmlands. The results showed that the NDVI spatial image information was complementary to, and as important as, the spectral information in models. In addition, biodiversity models generated using such reflectance-based data outperformed a model using the mean field size, the percentage of cropped area per landscape, and the the Shannon crop diversity index derived from a discrete land cover map that was generated by classification. Based on such empirical studies, and on the advantages of a reflectance-based analytical approach over the LULC map comparison for time-series analysis as described above, this study adopted a reflectance analysis approach.
Most remote sensing-based biodiversity studies have focused on the relationships of spatial variability in land cover or reflectance and biodiversity at a single point in time (e.g., [28,33,34]). However, in many ecosystems, and particularly in agricultural environments, temporal vegetation dynamics and their effects on farm species diversity may also be significant. Such vegetation dynamics may include changes in the amount of crop cover or density, as well as land cover changes such as the removal of natural vegetation in the non-cropped portions of the landscapes. The overall hypothesis of this study was that long-term inter-annual temporal trends and variability in farmland vegetation productivity over a large agricultural region are related to farm species diversity assessed at a given point in time.
More specifically, with respect to temporal vegetation trends, it was expected that vegetation productivity increased during the 2.5 decade study period, due primarily to agriculture intensification that has been noted in the region, with crops changing from less dense cereals to more dense corn and soybeans. However, the associated removal of small areas of natural vegetation at field margins that have occurred in some landscapes may temper such a trend.
With respect to the relationships of these vegetation dynamics with biodiversity, two aspects were investigated, vegetation productivity trends and variability. From the above general hypothesis, given significant positive trends in vegetation productivity were detected, the areal proportion of landscapes showing such trends was expected to be positively related to the diversity of all four taxa evaluated (see Section 2.1). More plant species were expected to be associated with increased vegetation productivity (NDVI), although the use of herbicides can temper that. Most of the beetle species prefer vegetation cover, so increased productivity was expected to benefit beetle diversity. We also expected birds to prefer more vegetation for nesting and food sources, and we expected pollinating butterflies to prefer more vegetation and plant species. With respect to the variability in long-term vegetation productivity, there is very little knowledge on which to base specific hypotheses. Only a generalized hypothesis could be made-that landscapes with higher variability due to intensive crop rotation, fields left bare or fallow in some years, and removal of natural vegetation at field margins in specific years, should have lower biodiversity in all four taxa due to difficultly for species to establish and maintain life cycle activities in such highly variable areas.
This study used time series of remote sensing data and statistical techniques to test the above hypotheses. The objectives of the study were to 1) determine if agricultural landscapes of eastern Ontario experienced inter-annual trends in vegetation productivity during the 2.5 decades prior to 2011-2012 (the years of biodiversity data collection), and 2) if such trends or temporal variability are related to 2011-2012 biodiversity, meaning such antecedent farmland vegetation dynamics potentially impact biodiversity. To address these objectives, a regional landscape-based approach was taken comparing similar spectral data at two scales of analysis using Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat data, each having rich temporal archives but with distinctly different spatial resolutions.

Study Site and Data
The study area is located within eastern Ontario, south of the City of Ottawa (Figure 1). It is approximately 10,000 km 2 and is dominated by agricultural lands that, as of 2011, were comprised mainly of corn (21%), soybeans (19%), forage crops (alfalfa, clover, hay; 30%), and wheat (3%) [35]. Prior to this study, 93 sample landscapes, as shown in Figure 1, were selected using a protocol [36] designed to maximize the variance in crop type composition and configuration heterogeneity amongst the sample landscapes. The number of agriculture land cover in selected sample landscapes was controlled to be between 60% and 90%, and the percentage agriculture in each sample landscape was found to not be related to the biodiversity indices described below (p > 0.05). Spatial autocorrelation between sample landscapes was also controlled by ensuring they were separated by a minimum distance that produced non-significant local Moran's I values [36]. Integrated analysis of the above factors was conducted for potential landscape sample sizes of 1 × 1 km to 7.5 × 7.5 km. The sample landscape extent that best satisfied the above criteria was 3 × 3 km. Biodiversity was surveyed prior to this study in the cropped portions of the sample landscapes (46 landscapes in 2011; 47 landscapes in 2012), i.e., species inventory was conducted within cropped fields. Logistically, biodiversity surveys could not be conducted across the full 3 × 3 km extent of each sample landscape. Instead, they were completed in the central 1 × 1 km extent of each sample landscape at four points that were at least 200 m apart and 50 m from both non-crop areas and the edges of the sample landscape. Seven taxa groups were surveyed, including birds, butterflies, bees, spiders, beetles, syrphids, and non-crop plants. Alpha diversity (the average number of species recorded at all sample locations within the landscape), beta diversity (the degree of changes in species composition and abundance between sample locations in the landscape), and gamma diversity (the total number of species over all the sample locations in the landscape and equal to the sum of alpha and beta) were calculated for each taxon. Details of the biodiversity methods and the list of surveyed species are given in Fahrig et al. [20]. The land cover within these 1 km × 1 km areas was also mapped in the field.
Prior to the research of this paper, initial exploratory statistical analysis ( [37]; not presented here for brevity), demonstrated that these diversity indices for birds, butterflies, beetles and non-crop (wild) plants showed some significant bivariate correlations (p < 0.05) with certain phenological variables and mid-summer NDVI derived from 2011 MODIS data as well as mid-summer 2011 Landsat 5 NDVI and Tasseled Cap components [38]. Based on this, these four taxa were retained for further analysis. In addition, alpha diversity of the selected taxa was not as strongly related to these remote sensing variables as were beta and gamma diversity. As correlations of gamma with alpha and beta were high (|r| > 0.88), and because gamma diversity is a simple indicator of overall species diversity at the landscape scale, it was selected for further analysis. The number of species surveyed over the 93 landscapes in each of the four taxa was as follows: 227 (plants), 80 (beetles), 52 (birds); and 30 (butterflies). The habitat preferences for all species within the types of farm fields/cover types where the surveys were conducted were not known in detail; typically, general information was available such as whether they preferred vegetation (most species), more open conditions, or if they were known to be generalist species. Thus, it was expected that the hypothesized positive vegetation trends in the sample landscapes would benefit most species in the four taxa. In addition, the scale of response to vegetation change was not known for all species, and therefore, only broad generalizations by taxa could be made as described previously in the study hypotheses.
The lack of detailed knowledge of the scale of response to vegetation change for all species also impacted the remote sensing data type selection. Higher resolution imagery would have been desirable to potentially detect smaller area vegetation changes to which species may respond, but such time series do not exist for the study area since high-resolution sensor data from Ikonos, Quickbird, and their successors are not acquired universally and archived. Such time series would also have been prohibitively expensive over the large study region. For these reasons, available time series archives of more than a decade at moderate (Landsat) and coarse (MODIS) resolutions were selected to compare their capability to detect temporal vegetation trends in the sample landscapes and model biodiversity. As a common and well-calibrated reflectance metric representing vegetation quantity or cover, mid-summer NDVI was selected to compare between these sensors and scales. As the overall hypothesis was that vegetation dynamics antecedent to the time of biodiversity measurement (2011-2012) are related to biodiversity, remote sensing data up to 2011 were selected (also, Landsat 5 data were not available after November 2011).
Time series data of MODIS weekly NDVI with 250 m spatial resolution were provided by the Canadian Ag-Land Monitoring System (CALMS) of the AAFC. CALMS is implemented throughout Canada's growing season (start of Julian week 12 (end of March) to end of Julian week 43 (late October)). It includes atmospheric correction and removal of clouds and cloud shadows in the creation of seven-day composites of the MOD09GQ Surface Reflectance (Davidson, 2015). In this study, the average of the July MODIS weekly NDVI from 2000 to 2011 was used, corresponding to the middle of the growing season.
Sixteen Landsat 5 Level-1 Surface Reflectance scenes (row 29/path 15) from 1985 to 2011 were available that had less than 5% cloud cover. Each had been atmospherically corrected by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS). Fourteen of the dates (Table 1) were within the mid-summer portion of the growing season (July to mid-August). Two scenes were outside this range by a small amount (August 20, 1987, and August 28, 1996) but still well before most crop harvesting and vegetation senescence, which typically occurs in mid-to late September. In addition,~20 cm pixel air photos were acquired over half the landscapes in 2011 and the other half in 2012. They were used to interpret and map land cover within the landscapes to aid validation of the results of the Landsat and MODIS time series trend analyses.

Linear Trend Analysis
Trend analysis was conducted to quantify inter-annual vegetation trends in the sample landscapes and investigate the potential impacts of those trends on biodiversity. The Theil-Sen (TS) slope linear trend measure [39,40] was selected to compute the trend magnitudes since it avoids spatial and temporal autocorrelation and is robust to outliers. It estimates the median of all n(n − 1) / 2 slopes calculated between observation values at all pair-wise time steps [41,42]. The Contextual Mann-Kendall (CMK) test [43] was used to assess the significance of the identified trends. The null hypothesis is that no trend exists, i.e., the data are randomly ordered and are not time-dependent. The principle of this method is to characterize ground features based on their spatial autocorrelation, presuming that a pixel would not show a drastically different trend from the neighboring pixels [42,44]. If an isolated location indicates a trend while its neighboring locations do not, it is expected that the trend is spurious. The CMK procedure tests the slopes between all pair-wise combinations of samples. The data are referenced to time, and each data point is the reference to the next data point in the sequential time period [43]. Furthermore, it directs all n(n-1)/2 iterations in a way that, if the value of a later data point is higher than an earlier one (i.e., the slope is positive), then the S statistic increases by 1, and if the opposite is true, it decreases by 1 [45]. The Kendall S is defined as: where n is the length of the time series, and xi and xj are the observations at times i and j, respectively. The variance of the S statistic, in relation to n is used to determine a significance (p) value, which indicates the probability of obtaining a value greater than or equal to S when no trend is present. Details of integrating the TS and CMK tests can be found in Neeti and Eastman [43]. The TS slope and CMK significance of trends were calculated for each MODIS pixel and each Landsat pixel within all 3 × 3 km sample landscapes. The percentage of pixels in each sample landscape that demonstrated a significant positive slope or a significant negative slope (p ≤ 0.05), was calculated for each 3 × 3 km sample landscape. The output of MODIS NDVI trend analysis (map of pixels with significant positive and negative trends) was resampled using the nearest neighbor approach to match the 30 m Landsat pixels and the maps were then cross-tabulated to calculate the numbers of coincident pixels of both datasets that demonstrated significant negative or positive trends, or no trend.

Biodiversity Modeling
To assess the relationships of vegetation trends and temporal variability with biodiversity, bi-variate Pearson linear correlation analysis was used. The gamma diversity of birds, beetles, butterflies, and wild plants was separately modeled against individual temporal variables from two groups. In the first group representing the overall trend over the given time series, models were created using the percentage of pixels in each sample landscape with a significant negative trend or positive trend, respectively. In the second group, predictor variables were based on the inter-annual variation in sample landscape vegetation as measured by the coefficient of variation (CV) of sample landscape NDVI over the time series. This temporal CV of NDVI was calculated for 1) the sample landscape spatial average of NDVI, and 2) the sample landscape spatial CV of NDVI. These two metrics represent the temporal variability in the mean sample landscape NDVI and the temporal variability of sample landscape spatial variability in NDVI. Bi-variate models created using the above variables derived from MODIS and Landsat were then compared. In addition, for taxa showing significant relationships with more than one explanatory variable, forward stepwise multiple regression (SWR) was implemented to determine if stronger models could be generated as well as the relative contributions of the variables to the models. The threshold for variables to enter the SWR model was p ≤ 0.05. Model residuals were checked for normality and uniformity, and none were found to violate these conditions.

MODIS and Landsat Midsummer NDVI Inter-Annual Trends
The evaluation of the TS slope and CMK significance for 2000-2011 July average MODIS NDVI data within the sample landscapes shows that 38% of all pixels demonstrated a statistically significant trend in NDVI. Many more pixels showed positive or negative trends that were not statistically significant [37]. Figure 2 shows the statistically significant positive (green) and negative (red) trends over the whole study region with the sample landscapes delineated. It is evident that there had been large areas of significant vegetation decline in the Ottawa region. This is generally due to urban development. These areas were not suitable for landscape selection because they were no longer completely agricultural in 2011-2012. Figure 3 displays close-ups of the results for three sample landscapes as an example. In Table 2, the percentage of pixels in each sample landscape that had significant negative or positive trends is shown as three groups (0%-4%, 5%-9%, and over 10%), along with the number of sample landscapes associated with each of these three groups. It is evident that more sample landscapes had greater proportions of pixels with significant positive trends in NDVI than negative trends. In particular, close to 1/3 of the sample landscapes (39 of 93) had more than 10% of their pixels showing a significant positive trend in NDVI over this period.
The results described above for MODIS NDVI data were more pronounced for the Landsat data, most likely because Landsat's temporal archive (1985-2011) was longer than for MODIS. A larger proportion of Landsat pixels (48%) showed a significant trend (Figure 4), and many more sample landscapes had larger proportions of positively trending pixels than for the MODIS data (e.g., more than 10% of pixels demonstrated significant positive trends in 42% of the sample landscapes for the MODIS dataset and in 96% of the sample landscapes for the Landsat dataset, Table 2). Almost all sample landscapes (90 of 93) showed significant positive trends in more than 10% of their pixels, and the highest percentage of pixels in a sample landscape with a significantly positive trend was 45%. The MODIS data were somewhat better at detecting larger area positive NDVI trends associated with crop fields, particularly in large fields. Overall, the Landsat data were found to be much better for the detection of vegetation trends within the 3 × 3 km sample landscapes of this study; MODIS data would be appropriate for more generalized detection of trends over larger spatial extents.

Detection of a Change in Trend Slope for Subperiods in the Landsat Time Series
As a consequence of the analysis of the spatial average of the Landsat NDVI data for each sample landscape over the full 1985-2011 period, it was observed that more than half the sample landscapes (56 of 93) demonstrated a break in the time series and a change in the direction of the slope of the trend line from positive in the earlier subperiod to a non-significant or sometimes significantly negative slope in the later subperiod. This break was most evident for these sample landscapes for the subperiods of 1985-2000 and 2000-2011, although in some sample landscapes, later years such as 2005 were stronger breakpoints. An example of one sample landscape (coded H_L_40) is presented in Figure 5. Given this observed change in trend slope, and that the year 2000 was also the start of the MODIS time series, these two time periods were selected to compare the early versus late trends in Landsat NDVI and to compare Landsat and MODIS trends for all sample landscapes (Section 3.2).  Comparing CMK significance for all pixels over all the 93 sample landscapes for 1985-2000 and 2000-2011 revealed that many more pixels displayed a positive trend in the earlier period compared to the later period. Figure 6 displays these results for both periods as five groupings, from 0-9% to over 40% of all pixels over all sample landscapes. Analysis at the landscape level showed that, for 1985-2000, the number of sample landscapes with a given percentage of pixels showing significant positive trends was as follows: 0%-19%-seven sample landscapes; 20%-29%-26 sample landscapes; 30%-39%-41 sample landscapes; >40%-19 sample landscapes. The highest percentage of significantly positive trending pixels in a sample landscape was 58%. In contrast, for 2000-2011, only 17 sample landscapes had more than 10% of their pixels, showing a significant positive trend (highest percentage for a sample landscape was 33.5%). Negative trending pixels were far fewer, being less than 5% in all the sample landscapes for both periods. These results indicate that most sample landscapes showed distinct positive trends in vegetation productivity in up to 58% of their pixels during the earlier period, while in the later period, these trends leveled off to become non-significant, or for a small number of pixels in the given sample landscapes, NDVI declined.  Figure 7 shows the results in map form for significant positive and negative trends over the whole study region and for some sample landscapes as an example. The total number of pixels over all the sample landscapes that displayed a significant negative slope was small in both periods but greater in the later period than in the earlier period. The negatively trending pixels in 2000-2011 were often associated with the removal of natural and semi-natural patches of trees or hedgerows (Figure 8). This was visually validated using of high resolution Google Earth imagery. In addition, most of the positive trends were associated with cropped lands, indicating increased green vegetation (Figure 9). This might be due to changes in agricultural practices or in landcover from crops with lower mid-summer NDVI to crops with higher mid-summer NDVI.

Comparing 2000-2011 MODIS and Landsat Mid-Summer NDVI Trends
The cross-tabulation of the re-sampled MODIS 2000-2011 trend maps with their corresponding Landsat trend maps found that 92% of the MODIS pixels showing no trend were in agreement with those of Landsat, and 88.4% of Landsat pixels with no trend were in agreement with those of MODIS. However, only 3.1% of the MODIS significant negative trending pixels and 6.3% of the MODIS significant positive trending pixels were in agreement with the coincident Landsat pixels, and only 6.3% of the Landsat significant negative pixels and 8% of Landsat significant positive pixels were in agreement with the coincident MODIS pixels (Table 3). Despite this, the general spatial distribution of these pixels coincided well in many areas. Landsat generally identified more pixels with significant positive or negative trends and provided more spatial detail. MODIS identified fewer areas of significant change, and the spatial extent of those change areas, was often over-estimated. Figure 10 shows two example areas with spatially concurrent significant trends in both datasets as ellipses, while areas with significant trends identified by only one dataset are indicated by arrows. This figure also shows that Landsat pixels with significant trends are more spatially precise since they mostly fall within the boundaries of a given vegetation unit (field, forest, etc.). The coarse MODIS pixels often extended outside such boundaries in areas where change had not occurred, and they often extended outside the given sample landscape.  Overall, the Landsat data were found to be much better for the detection of vegetation trends within the 3 × 3 km sample landscapes. MODIS data were somewhat better at detecting larger area positive NDVI trends associated with crop fields, particularly in large fields, and such data would be appropriate for more generalized detection of vegetation trends over larger spatial extents than the 3 × 3 km landscapes of this study.  (Figure 11). In Figure 11a (left), there are three distinct sample landscapes separate from the other data points at the upper end of the models; these sample landscapes had more than 20% of their pixels with significant positive NDVI trends. Removal of these three points reduced the correlation, but it was still significant (r = 0.33, p = 0.00); hence, these sample landscapes were not deemed to have a dominating effect on the relationship found.
Overall, these results show that where there was a gradual but significant increase in vegetation productivity during 2000-2011, plant species diversity measured over the sample landscapes at the end of this period was higher. This indicates that it is the decade preceding the biodiversity measurements that is more important than the decade before that in influencing biodiversity.

Relationships between Biodiversity and Inter-Annual Variability in Vegetation Productivity
In biodiversity modeling using the inter-annual variability of NDVI at the landscape level (as represented by the temporal CV of the mean sample landscape NDVI and of the spatial CV of each sample landscape), no significant relationships were found for the MODIS data. Significant negative relationships (p < 0.01) were found between the 2011-2012 plant and beetle gamma diversity and Landsat temporal CV of the mean sample landscape NDVI for all the three-time periods as follows (Figure 11b

Multiple Regression of Plant Gamma Diversity vs. Landscape Temporal Dynamics Metrics
Plants were the only taxon that demonstrated a significant relationship of gamma diversity with both types of temporal variables (temporal CV of the mean sample landscape NDVI; the percentage of significant positively trending pixels in each sample landscape), and this occurred for the 2000-2011 Landsat NDVI time series. The plant gamma diversity of 2011-2012 was, therefore, used in an SWR model against these two 2000-2011 temporal variables to test the relative importance of each variable in the model. These variables were not correlated with each other (r = −0.31) so multicollinearity was not an issue.
The resulting two-variable model demonstrated a stronger relationship (Adjusted R 2 = 0.25, p = 0.00) than the best single variable model R 2 for the 2000-2011 period. The percentage of significant positive trended NDVI pixels was the first variable entered (r = 0.44), while the temporal CV of the mean sample landscape NDVI was the second variable entered (r = −0.27). The results show that both temporal variable types are complementary, and, together, they strengthen the evidence of the initial hypothesis that temporal vegetation dynamics prior to a given year when biodiversity is assessed are significantly related to and potentially have an impact on biodiversity.

Discussion
This study had two main components: 1) detect and evaluate temporal trends in vegetation productivity in agricultural landscapes of eastern Ontario during the years preceding 2012 using MODIS and Landsat data, and 2) determine if these vegetation dynamics are related to farmland biodiversity measured in 2011-2012.

Trends in Vegetation Productivity
Significant positive trends in the mid-summer NDVI data were found in all sample landscapes with the proportion of pixels showing such trends being over 1/3 for MODIS (2000-2011) and almost 1 /2 for Landsat . Most of these trends occurred in the cropped portion of the sample landscapes, indicating increased green vegetation due either to agricultural practices or changes from crop types, such as alfalfa and winter wheat, which are harvested in this region in summer, to denser foliage crops such as corn and soybean; demand for the latter two crops increased in the last two decades for feed (both) and ethanol fuel (corn) [46]. According to the Census of Agriculture produced by Statistics Canada [2], a major shift in crop composition in Ontario in the early 2000s was reported from less heavily managed crops such as hay, alfalfa, and perennial plants to more intensively managed crops such as corn and soybean. These positive trends were most evident in low heterogeneity sample landscapes with few crop types and large fields [37].
While the percentage of pixels showing significant negative trends in sample landscapes was smaller, those pixels were often spatially clustered and associated with the removal of woody vegetation or destruction of natural patches for either crop production or new buildings. Higher proportions of pixels in sample landscapes with high compositional heterogeneity (more crop types) and configurational heterogeneity (smaller field sizes), both of which are known to benefit biodiversity [20], showed negative trends in vegetation productivity [37] in a general process of agricultural intensification over the study period. The detection of smaller negative trending areas within sample landscapes was facilitated much better using the higher spatial resolution Landsat data. Although up to 144 MODIS pixels overlapped with each sample landscape (depending on pixel alignment and image orientation with respect to sample landscape orientation), at 250 m pixel spacing, they were too large to consistently detect removal of trees and hedgerows. This concurs with Rogers et al. [47] who found coarser spatial resolution AVHRR and MODIS time series to be poorer than Landsat time series in portraying temporal changes and detecting early signs of tree mortality in the boreal forests of North America.
Additionally, during 1985-2000, many Landsat pixels in the sample landscapes demonstrated significant positive trends while between 2000 and 2011, many pixels demonstrated no trend or a significant negative trend. The analysis of such sub-periods could be extended to shorter intervals (e.g., 2-4 years) to evaluate the effects of crop rotation, other landscape changes, or short-term weather variations. Different time intervals could be selected to test the effect of time series length on the trends, to identify different sub-trends, and to test the sensitivity of trends to the selected time intervals between images. For example, Czerwinski et al. [44] divided a 20-year time series of Landsat 5 images into 1-3, 4-5, and 6-7 year intervals to analyze time series extent effects on the detection of forest changes, finding that the 1-3 and 4-5 year intervals produced similar change results, while the 6-7 year interval did not detect the same changes and suggested very little changes had occurred.

Relationships of Biodiversity with Inter-Annual Landscape Vegetation Dynamics
As for the trend analysis, Landsat data were much better for modeling biodiversity than MODIS data (which did not produce any significant models). For Landsat data, plant gamma diversity measured in 2011-2012 was found to be negatively related to the temporal variability (CV) in NDVI during the preceding 26 years, as well as for the two sub-periods of 1985-2000 and 2000-2011. Plant diversity was also positively related to the proportion of pixels with a significant positive trend for 2000-2011 and 1985-2011 in each sample landscape. However, this relationship was stronger in 2000-2011 compared to 1985-2011, indicating that a more recent trend in vegetation may have a stronger influence on plant diversity. In the stepwise regression, a two-variable model for plant diversity was derived with the percentage of pixels with a significant positive NDVI trend accounting for more variance than the inter-annual variation (CV) in NDVI. These results for plants were generally as expected. In sample landscapes with more pixels showing a long-term increase in vegetation quantity and where inter-annual variability in vegetation quantity is lower (i.e., less year-to-year variability in NDVI), more wild plant species could be established. For beetles, higher inter-annual variation in vegetation quantity was also associated with reduced gamma diversity. Most of the beetle species surveyed prefer vegetated habitats. Therefore, high inter-annual variation in the mid-summer vegetation quantity, where in some years vegetation is thinner with greater soil exposure due to weather, the selected crop type, and phenology, or other management practices, results in lower beetle diversity.
Despite being statistically significant, these relationships account for low amounts of variance in gamma diversity, i.e., they are noisy. This may be because a process such as agricultural intensification that has occurred within the study area during the decades of the data used for this study may have contrasting effects. The relationships found indicated that a greater proportion of pixels with positive long-term NDVI trends in a sample landscape is associated with higher plant diversity, and this trend in sample landscapes is most likely due to transitions to higher biomass crops such as corn and soybeans. However, through such intensification, many farms have also experienced the removal of natural vegetation in woodlots and at field boundaries, which reduces landscape NDVI, while herbicide use increases, directly killing off some non-crop plant species. The former increases noise in the x-variable (% pixels with significant positive trend) while the latter increases noise in the y-variable (plant gamma diversity). The relationships found here are, therefore, not appropriate for the prediction of biodiversity, but they do give statistical evidence that plant and beetle diversity are related to vegetation trends over the long-term as well as inter-annual vegetation dynamics. This is important new knowledge with respect to the development of an understanding of impacts of farming on biodiversity in the temporal dimension; it is complementary to the multitude of studies that have shown impacts of landscape spatial heterogeneity on biodiversity. This is one of the first, or the first study evaluating effects of such temporal trends and inter-annual vegetation dynamics on biodiversity in agricultural landscapes using remote sensing time-series data. Previous studies found negative impacts of agricultural land use change. For instance, Senapathi et al. [23] used digitized land cover maps to assess the impacts of 80 years of land cover change on bee and wasp species richness and community composition in England. Reidsma et al. [22] used the European Union Farm Accountancy Data Network to classify different farm types and land use intensity, and assessed and projected the impacts of agricultural land use change on biodiversity for 2010, 2020, and 2030 scenarios. Jiang et al. [21] used data from the ground sampling of three sites with different land use practices to assess the impact on plant species richness. However, remote sensing-based time series data were not utilized in those studies and, therefore, studies in other environments and contexts were consulted for comparison. At broader scales, De Baan et al. [48] found generally negative impacts of change in different land use types on biodiversity in various regions of the earth. Jetz et al. [49] estimated the projected impacts of climate change and land use change on global avian diversity and found that 400 species are projected to suffer by the year 2050. Such studies have a temporal dimension but generally use land cover maps created for one or more years instead of a time series of remotely sensed reflectance data, as in this study.

Limitations and Recommendations for Future Research
This study was constrained by the biodiversity data, which were acquired in 2011-2012 as part of a large field and remote sensing project (e.g., [20,33,36] amongst others). The acquisition of the biodiversity data was a major undertaking by teams of up to 10 field personnel over the April-September period in both years (1/2 the sample landscapes surveyed in each year). The use of such intensive and directly acquired biodiversity data provides an excellent reference for such remote sensing modeling studies. However, being essentially at one point in time limits such a temporal study as this to the analysis of the potential effects of antecedent vegetation dynamics on biodiversity measured at that time. While the results of this study are illuminating and point to the potential impacts of temporal landscape dynamics on farmland biodiversity, better temporal biodiversity data could be pursued, using continuous multi-year field efforts. In fact, a portion of the sample landscapes of this study were revisited in 2016, and 2017 and bee and butterfly species were surveyed using the same protocols by a field crew in the Landscape Science and Technology Division in Environment and Climate Change Canada [50]. The sample landscapes may be visited in the future and used as long-term monitoring sites to assess the effects of changes in the landscape as well as climate factors on pollinator diversity in Eastern Ontario. In addition, other biodiversity indicators may be acquired over a given time period. However, Hanley et al. [4], in a study using density-based metrics and species richness in models based on farm types with various land uses (e.g., grazing and sheep, dairy, and beef production), management practices and policy scenarios, caution on using umbrella species or generalized biodiversity indicators such as total richness, which can miss important species compositional effects. Sources such as the long-term USGS North American Breeding Bird Survey [51] are recommended for use as a temporal reference dataset for modeling against landscape dynamics as derived from remote sensing time series. In addition, comparing farmland biodiversity with the biodiversity of adjacent natural and semi-natural lands is recommended to assess the changes in species diversity due to land conversion to agriculture.
With respect to the remote sensing data used in this study, the coarse spatial resolution of the MODIS data was expected to be a hindrance in trend analysis and biodiversity modeling in comparison to Landsat data for the given landscape scale. However, its excellent temporal resolution allowed better temporal precision of average mid-summer NDVI estimation than for Landsat, where single images had to be used that varied in acquisition date, potentially adding noise to the trends and biodiversity models. Thus, MODIS is more useful for larger regional or global inter-and intra-annual trend analysis, as others have found. Alternatively, fusion techniques could be investigated to take advantage of the MODIS high temporal resolution and Landsat high spatial resolution (e.g., [52][53][54] in studies of landscape changes and trends). High-resolution image time-series data from other sensor platforms (e.g., air photos; Quickbird-Worldview) were not archived for the study area and would have been very costly to acquire over the spatial and temporal extents of this study. Where the resources exist to acquire such data over both a suitably long period and with sufficient areal coverage to represent a given eco-zone (such as the farmland region of this study), future research should consider such data. They should be able to detect smaller areas of vegetation trends with better spatial precision than Landsat, and some taxa may respond to such small vegetation changes. However, many or most species in the four taxa of this study are known to function and disperse over larger areas (e.g., by flight or wind), so it is uncertain if higher-resolution detection of vegetation trends and variability in small areas would improve such biodiversity modeling. Such data could be used to try to determine the scale of the effects of vegetation dynamics on biodiversity to see if there are differences among taxa. This is an interesting question that we could not tackle with our datasets. In addition, comparison of LULC classifications derived from such high-resolution imagery over time could aid in the identification of the causes of the detected vegetation trends and variability, provided error propagation is not significant as discussed in the Introduction. Analysis of patch-based metrics derived from such maps may help elucidate processes such as changes in connectivity over time that could affect the biodiversity of certain species or taxa. These could be useful in guiding the development of policy and guidelines regarding farmland management over time for the enhancement of the diversity of given taxa.

Implications of This Study for Policy and Monitoring of Biodiversity
In the past few decades, Canadian agricultural policy has shifted from increasing output and promotion of income stability in the agriculture sector to research focusing on sustainable ways to increase productivity and limit environmental impacts [55]. To fulfill this, Agriculture and Agri-food Canada (AAFC) uses a series of agri-environmental indicators in the analysis of overall environmental performance. These indicators are assessed over time, and areas showing improvements are highlighted, while regions that require policy interventions are identified. For example, the Wildlife Habitat on Farmland Index [55], as an indicator of biodiversity, demonstrated a decline from a "moderate' state in 1986 to a "poor" state in 2011 in Canada.
Indices such as this that are developed and employed at the national or provincial level are often based on mapping of habitat using field observations and classification of moderate to coarse resolution remote sensing imagery. Given that comparison of classifications from different years and over a time series is sensitive to error propagation as mentioned previously, aggregated amounts of 'habitat' are often generated by year and compared between years, as opposed to specific statistics on class trends and changes from one class to another. The main advantage of the methods of this study is that they do not require such a classification step and avoid the associated map error propagation in the time series analysis. The methods of this study could be automated and used operationally to evaluate vegetation productivity trends and variability as a type of index of potential biodiversity (at least for plants and beetles, as shown here). In addition, a hybrid of the classification and reflectance approaches could also be developed, for example, by delineating areas of significant NDVI trends or high variability using segmentation. Indices of the abundance and spatial properties (e.g., fragmentation, connectivity metrics) of these identified patches could then be derived and used in more detailed landscape vegetation and habitat analysis that feed into policy interventions of the type mentioned above for AAFC. The demonstration and justification to farmers of the value of maintaining stable crop field vegetation cover from year to year and the preservation or installation of natural and seminatural vegetation are suggested as means to enhance farmland biodiversity and reinforce such agricultural policies.

Conclusions
This study evaluated relationships between long-term farmland temporal vegetation productivity dynamics and 2011-2012 biodiversity of four taxa (plants, beetles, birds, and butterflies). It compared the potential of MODIS and Landsat 2000-2011 time-series datasets, as well as Landsat for the additional periods of 1985-2000 and 1985-2011, in detecting significant positive or negative trends in NDVI in 93 sample landscapes of 3 x 3 km in a predominantly agricultural area of eastern Ontario. Over all sample landscapes, 38% of MODIS (2000-2011) pixels, and 48% of Landsat (1985-2011) pixels showed significant positive or negative NDVI trends, with positive trends being more dominant than negative trends. Most of the significant negative trends were spatially associated with removal of natural and semi-natural patches, and, in some cases, the removal of field edges and hedgerows, and these were better detected with the higher resolution Landsat data. The most significant positive trends in NDVI were in crop fields and probably associated with increasing crop density, either by changes in crop type over time or by management practices to increase growth and yield. The full 1985-2011 period of Landsat data showed a break in the slope of NDVI about the year 2000 for many sample landscapes where a positive slope from 1985-2000 changed to a non-significant or, in some cases, a negative slope for 2000-2011.
While MODIS data did not produce significant models of biodiversity, using Landsat data, the proportion of sample landscape pixels with significant positive trends and the inter-annual variation (CV) of the mean sample landscape NDVI were significantly related to the 2011-2012 gamma diversity of plants (positively and negatively, respectively). For beetles, gamma diversity was also significantly and negatively related to the temporal CV of the mean sample landscape NDVI. The relationships for the CV of NDVI were stronger for the overall Landsat data period (1985-2011) than for either of the sub-periods of 1985-2000 and 2000-2011. In future work, higher resolution airborne or satellite image time series, if available and affordable for a long-term regional study such as this, would help identify smaller areas of vegetation dynamics that may be more strongly associated with species and taxa life cycle functions, thereby potentially providing stronger biodiversity models.
Besides this study, the effects of inter-annual changes in remote sensing-based agricultural landscape vegetation indices on biodiversity have not been largely studied. The results found here show that long-term trends and inter-annual vegetation variability are complementary indicators of the potential effects on farmland biodiversity and can be used to inform landscape management policy and practice. This study has also demonstrated the value of remotely sensed reflectance time series data in the detection of vegetation trends and the modeling of biodiversity. This approach is a viable and simpler alternative to the classification and temporal thematic map comparison analysis. It depends only on well-calibrated reflectance data time series, which are freely available for several sensors, in contrast to temporal analysis of thematic classifications, which is highly sensitive to error propagation. Further research is warranted using other temporal analysis techniques and additional temporal metrics as well as temporal biodiversity datasets.