Arthropod biomass increase in spring correlates with NDVI in grassland habitat

Data from remote sensing are often used as proxies to quantify biological processes, especially at large geographical scales. The normalized difference vegetation index (NDVI) is the most frequently used proxy for primary productivity. Assuming a direct, positive interrelation between primary and secondary production in terrestrial habitats, NDVI is often used to predict food availability for higher trophic levels. However, the relationship between NDVI and arthropod biomass has rarely been tested. In this study, we analyzed extensive datasets of arthropod communities from semi-natural grasslands in central Europe to test the relationship between arthropod biomass of consumer trophic levels (“herbivores,” “mixed,” and “carnivores”) in grassland communities and NDVI during the spring season. We found that arthropod biomass generally increased with NDVI. The same positive relationship between biomass and NDVI was observed for each individual trophic level. Cross-correlation analyses did not show statistically significant lags between the NDVI and biomass of herbivores and carnivores. All in all, our study provides correlational evidence for the positive relation of primary and secondary productivity in temperate terrestrial habitats during spring. Moreover, it supports the applicability of NDVI data as a suitable habitat-specific proxy for the food availability of insectivores during spring.


Introduction
The availability of nutritional resources can critically influence community composition and species richness (Mittelbach et al. 2001;Ribas et al. 2003). Animals in seasonal habitats must deal with temporal variation in resource availabilities caused by periodic changes in vegetation structure (Stinson and Brown 1983) and habitat quality (Wiegand et al. 2008). However, estimating resource availability is often difficult to achieve, especially when large and/or remote areas should be sampled. Thus, ecologists frequently use proxies derived from remotely sensed data (Pettorelli et al. 2011;Stephens et al. 2015). Their application requires knowledge about the underlying relationship between the measured and the predicted traits. Such interrelation is relatively easy to establish if these traits are directly linked, as occurs between marine primary productivity and phytoplankton abundance (Marshall and Nesius 1996) or primary productivity and plant biomass in terrestrial ecosystems (Scurlock et al. 2002). However, approximations encompassing multiple trophic levels are harder to evaluate.
In ecology, a frequently used indicator for resource approximations is the normalized difference vegetation index (NDVI), a measure of primary productivity used for many years in remote sensing. NDVI distinguishes the reflectance properties of the vegetation (NIR-Red)/(NIR + Red), with NIR representing the near-infrared and red the visible red light reflected by the surface (Myneni et al. 1995). Nowadays, NDVI data are globally available (e.g., www.star.nesdis. noaa.gov). The temporal and spatial resolutions of these data have improved over time and are, therefore, ideal for the analysis of vegetation dynamics at large scales.
The applicability of NDVI for quantifying various habitat properties has been repeatedly demonstrated. Thus, NDVI has been proven as a good proxy to estimate plant biomass (Weiser et al. 1984), vegetation quality (Hamel et al. 2009), and crop yield (Fuller 1998). Moreover, NDVI has been used to predict food availability for animals at higher trophic levels such as raptors overwintering in the Sahelian zone feeding on grasshoppers (Trierweiler et al. 2013) or kestrels at their breeding sites feeding on small mammals (Smith et al. 2017). Various predators like insectivorous mammals and birds adjust their annual cycles, especially the timing of breeding and migration (Perrins 1970;Jonzén et al. 2007), to the abundances of terrestrial arthropods, which serve as their main food source (Bell 1990;Morse 1971).
Several aspects concerning NDVI and arthropod biology have been revealed by previous research. For instance, Lassau and Hochuli (2008) demonstrated that NDVI of the open forest canopy is positively correlated with species richness and abundance of beetles in Australian dry forests. Changes in NDVI of birch forest has been shown to correlate with caterpillar abundance of geometrid moths in Fennoscandia (Jepsen et al. 2009), and NDVI of shrubs has been found to correlate with locust abundance (Deveson 2013). It is also known that the response of arthropod abundances to changes in NDVI can be delayed as a result of developmental time from larva to imago or reproduction of the next generation (Apiwathnasorn et al. 2006). Thus, the relation between NDVI and arthropods can differ according to life-history traits of species.
Arthropod biomass is a good parameter to quantify resources available for consumers of higher trophic levels. However, studies directly relating NDVI to arthropod biomass are scarce (Sweet et al. 2015). Thus, ecologists urgently need empirical support for using NDVI as a proxy for resource availability in insectivorous species.
Here, we related vegetation productivity derived from NDVI data to the biomass of arthropods in semi-natural grasslands in central Europe. Specifically, we investigated if and how the arthropod biomass of different trophic levels responds to the seasonal increase in NDVI that occurs during spring. As the phenology of arthropods may vary due to species-specific developmental times (Jarošík et al. 2011), we expect trophic levels to differ in their response to changes in vegetation phenology, i.e., the existence of lags between NDVI and biomass within the trophic levels.

Study area
We studied the relationship between primary productivity and arthropod biomass in semi-natural grassland (semi-dry meadows) at three sites (Gleitz, Johannisberg, and Leutratal) in eastern Germany between 1987 and 1989 (Appendix A).
The study sites (180-300 m a.s.l.) were situated on hill sites in distances of about 30 km. The regional climate is humidcontinental with low annual precipitation (< 600 mm p.a.) and late winter as the driest period during the year. During the 1980s, the thermic vegetation period (i.e., air temperatures above 5°C) averaged at about 230 days (Bundesamt für Naturschutz 2014).

Arthropod data
Arthropods were collected using identical methods throughout all three study sites and study periods within the framework of the grassland project run by the Institute of Ecology, Friedrich-Schiller University Jena (Germany).
In each study site, the sampling was carried out in two locations all year round, every year, using pitfall traps and sweep-nets for collecting ground and herb layer dwelling arthropods, respectively (for details see : Perner 1997;Malt and Perner 2002;Perner and Malt 2002). Within every location and every 2 weeks, ten pitfall traps were checked and ten double swings of sweep-netting in each of ten points were sampled. Each of these ten sampling points was randomly distributed coordinates around every location, in its vicinity, and with similar characteristics. Pitfall traps operated continuously, and the individuals trapped therein were collected the mid-day on the same day as sweep-netting was carried out.
In most cases, arthropods were determined to species level and, in some cases, only to family-level.

Trophic levels and biomass
We assigned each species to one of the following trophic levels: "herbivores," "mixed," and "carnivorous arthropods" based on information on main diet and foraging behavior (Table 1). "Mixed" category included omnivorous species and those that forage as herbivorous or carnivores during different stages of their life cycle. If species' specific information was not available (64% of cases), we categorized the species according to the typical trophic level of its genus or family.
To estimate the biomass per species, we multiplied species abundance (i.e., the total number of individuals irrespectively of the life stage) with body mass according to the following allometric relationship between body mass and body length: body mass (mg) = 0.0305 × body length (mm)^2.62 (Rogers et al. 1976). Data on typical body length for each species were obtained from Gossner et al. (2015). For the species not described in Gossner et al. (2015)-32% of the total-body length information was obtained from Watson and Dallwitz (2015) instead. Afterward, we summed biomasses of species within the trophic levels, for each sampling date and site, irrespectively of the method of collection. Finally, biomass data were log-transformed to achieve normal distribution (log (n + 1) to avoid negative values).

NDVI data
We acquired NDVI data from the weekly 16 × 16 km-gridded global dataset NESDIS STAR VHP 16 available by of the National Oceanic and Atmospheric Administration (NOAA) (https://www.star.nesdis.noaa.gov/smcd/emb/vci/VH/vh_ftp. php) using the R package rndvi from GitHub (https://github. com/tavimalara/rndvi). For each study site and sampling date, we assigned the NDVI value for the grid cell that includes the study site and corresponds to the closest date available in the NOAA dataset. Due to the close proximity between the study sites Johannisberg and Leutratal, the NDVI grid for the two sites was the same. We only considered the main growing season of vegetation, i.e., dates between the minimal and maximal values of NDVI for each year. Therefore, we excluded all NDVI and arthropod data after the annual NDVI maximum. The median duration for this interval was 12 weeks (25-75% percentiles: 11-14 weeks). For descriptive comparison between years, we determined the spring green-up date and the speed of green-up ( Fig. 1; Appendix E), which refers to the date for the maximum slope of NDVI rising in a certain year. These measurements inform about early vs. late spring seasons, and the speed of the primary productivity increase. Both measurements were calculated by fitting a logistic curve to the NDVI rise in spring, determined by its observed min and max values, for each site and year.

Statistics
Because arthropod biomass may respond with delay to changes in vegetation (Apiwathnasorn et al. 2006), we checked for time lags between NDVI data and the response in arthropod biomass using cross-correlation analysis of weekly mean NDVI and mean biomass. Time lags were considered relevant if the cross-correlation coefficient for each period was above a certain significance level. To improve the cross-correlation analysis, possible temporal autocorrelation within each time series were removed by using the package TSA (Chan and Ripley 2012).
We described the variability of NDVI and biomass data across sites, trophic levels, and years using LOESS curves. Further, we tested for the effect of NDVI-the explanatory variable-on arthropod biomass-the response variable-by applying a generalized linear model using the function lm of the package stats (R-Core Team 2018). We included trophic level as an explanatory variable and its interaction with NDVI (log(biomass)~NDVI + trophic level + NDVI * trophic level). Homoscedasticity and normal distribution of residuals were visually checked in the diagnostic plots. Fig. 1 Boxplot comparison (middle line shows the median, rectangles show 25-75% quantiles, and bar limits the range) of the two main studied variables in three study sites in Germany. a NDVI values at the study sites with different NDVI grid cells. b Arthropod biomass (g/sampling day) of species considered within "mixed" (including omnivorous species and those that forage as herbivorous or carnivores during different stages of their life cycle), "herbivores," and "carnivores" trophic levels collected at these sites We estimated the model coefficients using the maximum likelihood approach. All analyses were done in R (R-Core Team 2018).
We derived similar results when each trophic level had been analyzed separately. Herbivores' biomass was positively related to NDVI (t = 10.49; P < 0.001), carnivores' biomass also showed a positive relationship to NDVI (t = 5.34; P < 0.001), and biomass of "mixed" trophic level was positively related to NDVI as well (t = 2.72; P < 0.01). Assumptions for linear modeling were visually checked, and its requirements fulfilled.
Cross-correlation analysis did not show statistically significant cross-correlation coefficients for any lag considered between the NDVI and biomass of herbivores and carnivores. The lack of a lag suggests that biomass is related to the NDVI value measured most closely to the arthropod sampling date, as available in NOAA's weekly NDVI dataset.

Discussion
We found a direct, positive relationship between NDVI and arthropod biomass during the main growing season (i.e., period from the minimal to the maximal NDVI value of each year) in grassland habitats of eastern Germany. Thus, we provide (correlational) evidence for the interaction between a primary production index and the concomitant change of secondary production, which is frequently expected but has been rarely tested. The use of NDVI as a proxy for arthropod biomass allows for temporal assessments of secondary productivity in Table 2 Model output for the generalized linear model (GLM) relating arthropod biomass of the trophic levels ("mixed," "carnivores," and "herbivores") to NDVI, including its interaction. "Mixed" trophic level was used as reference level. The table shows the estimates of the resulting coefficients, their standard errors (SE), and test statistics to assess the rejection of the null hypothesis (under which these estimates are zero). Measures of goodness of fit (D-squared proportion, null deviance, and residual deviance) of the models are included as footnotes n.s. not significant Fig. 2 Relationship between arthropod biomass (g) (log-transformed) and NDVI. Points show observed summed biomass (g) for each site, date, and trophic level; the black line indicates the predicted relation obtained for the generalized linear model (log(biomass)~NDVI + trophic level + NDVI * trophic level), with 95% confidence intervals shown in grey shade various fields, e.g., phenology of food availability for insectivorous animals. We are confident that our conclusions are robust because our study covers a significant part of the local arthropod community, which had been collected in different vegetation strata by different methods as well as across various sites and years. However, our approach has some constraints which should be considered for generalization of the findings.
Our study focused on the growing season starting in spring, i.e., the spring green-up when seeds germinate, and new shoots and leaves grow (Cole et al. 2015), and when arthropods become active, i.e., after hatch or cease of the winter dormancy period due to rising ambient temperature (Jarošík et al. 2011). The beginning of the growing season offers food with the highest protein contents and digestibility to herbivorous arthropods (Crawley 1983;Mattson 1980), so it seems plausible that herbivores track this seasonal food availability by increasing abundances (Kerby et al. 2012). More active and abundant arthropods represent higher food availability for predators, which fine-tune the timing of arrival on the breeding grounds (Jonzén et al. 2007) and the start of reproduction according to spring phenology (Perrins 1970). Our study clearly indicates that NDVI is a reasonable predictor for arthropod abundance in temperate grassland habitats during the beginning of the growing season, and thus can be used as a proxy for food availability for organisms at higher trophic levels, e.g., migratory insectivorous birds at arrival or at the beginning of their breeding (Newton 2004;Grande et al. 2009).
However, our short study period of a few subsequent years prevents further conclusions regarding differences in maximum NDVI values between years and the related variation in arthropod biomass. Our study system is situated in the temperate climate zone where between-year variation in peak values is usually small and temporal shifts are more prominent (early vs. late years; see also Appendix E for differences between years). Moreover, NDVI data for our study period had a rather coarse resolution (16 × 16 km), and thus, a pixel representative for the study site also includes surrounding habitats, e.g., deciduous and coniferous forests, orchards, broad-leaved trees, or other habitats (Appendix B). Thus, the NDVI value represents a blend of habitat around the center of the satellite grid cell. We believe that this limitation is not likely to affect the analysis of temporal changes in NDVI, i.e., primary productivity in spring. Nevertheless, this limitation might compromise further comparisons of NDVI between sites and habitats during certain seasons, especially if surrounding habitats strongly differ between sites. If, for instance, forests remain green and photosynthetically active, NDVI values generally remain constant, while herbaceous vegetation might be already dried out, lowering the NDVI value of that cell. Moreover, annual maximum NDVI values of a certain habitat might not reflect maximum total primary productivity (Pinter et al. 1985;Wang et al. 2005) as various strata can differ in their seasonal pattern productivity, and certain groups of arthropods may still respond specifically to strata-specific variation.
In this study, we assumed that arthropods respond similarly to environmental conditions in spring. Moreover, we accounted for the vast variation in numbers of individuals and body masses of various species by categorizing distinct trophic levels. This procedure has a drawback for the applicability of NDVI-arthropod biomass interaction: to predict biomass quantities, researchers should account for differential trophic level abundances of arthropods species. However, the segregation of trophic levels allowed for checking for temporal lags in the interaction of biomass to NDVI rise. Such response lags are common as arthropods needs a certain time to complete their development; for instance, mosquito abundance rises 1 month after the first rainfall in the tropics (Apiwathnasorn et al. 2006). For temperate grasslands, we found no lagged response between rising NDVI and arthropod biomass of the two main trophic levels studied ("carnivores" and "herbivores"). Therefore, arthropods at our study sites more probably responded to NDVI value within the measurement time intervals of 2 weeks between subsequent NDVI measures. Thus, we contend that by using NDVI to predict arthropod biomass in temperate regions, developmental times must not be necessarily considered. However, trophic levels should preferably be considered when quantifying biomass from NDVI, as the shape of the positive relation between biomass and NDVI can vary for the different levels.

Conclusions
The positive and direct relationship between primary production and arthropod biomass confirmed by this study allows predictions of arthropod biomasses in spring to be made from remote sensing data, at least for temperate herbaceous habitats. Such assessments of secondary productivity can help to disentangle phenology and distribution patterns of consumers, i.e., the timing and the choice of stopover sites for birds on migration (Buler et al. 2007;Hutto 1985), the arrival and timing of breeding in insectivorous animals (Both et al. 2009), or may assist pest control management in agriculture at the beginning of the growing season (Bravo et al. 2003).
To generalize the relationship between NDVI and secondary productivity for various habitats, years, and locations, it is necessary to use long-term data sets of arthropods and detailed entomological descriptions. Additionally, the currently and freely available remote sensing data with higher resolution (e.g., 10 × 10 m resolution of Sentinel-2/ EU-Copernicus-Project) will enable future studies on primary productivity and related arthropod abundances across years, different management schemes, and on small scales.