Use of Carabids for the Post-Market Environmental Monitoring of Genetically Modified Crops

Post-market environmental monitoring (PMEM) of genetically modified (GM) crops is required by EU legislation and has been a subject of debate for many years; however, no consensus on the methodology to be used has been reached. We explored the suitability of carabid beetles as surrogates for the detection of unintended effects of GM crops in general PMEM surveillance. Our study combines data on carabid communities from five maize field trials in Central Europe. Altogether, 86 species and 58,304 individuals were collected. Modeling based on the gradual elimination of the least abundant species, or of the fewest categories of functional traits, showed that a trait-based analysis of the most common species may be suitable for PMEM. Species represented by fewer than 230 individuals (all localities combined) should be excluded and species with an abundance higher than 600 should be preserved for statistical analyses. Sixteen species, representing 15 categories of functional traits fulfill these criteria, are typical dominant inhabitants of agroecocoenoses in Central Europe, are easy to determine, and their functional classification is well known. The effect of sampling year is negligible when at least four samples are collected during maize development beginning from 1 April. The recommended methodology fulfills PMEM requirements, including applicability to large-scale use. However, suggested thresholds of carabid comparability should be verified before definitive conclusions are drawn.


Introduction
Although genetically modified (GM) crops are generally considered safe for non-target arthropods [1], there are still uncertainties regarding the long-term effects caused by the accumulation of miniscule changes in the agroecosystem. Science-based post-market environmental monitoring (PMEM) is therefore required for GM crops by EU legislation (Directive 2001/18/EC, [2,3]). PMEM aims to identify risks that did not become evident during the pre-market risk assessment. Consequently, PMEM results are expected to provide a basis for subsequent regulatory decisions, including the prolongation and modification of the monitoring plans. The detection of adverse changes in the environment may trigger additional research that could eventually lead to the withdrawal of approval for GM crops [4].
There has been much discussion about the PMEM of GM crops, but a general PMEM plan accepted by regulators, scientists, and the agricultural biotech industry is still lacking [5]. PMEMs are a legal

Characterization and Quantitative Comparison of Carabid Communities
The sum of catches in all localities totaled 58,304 individuals belonging to 86 species. Within functional traits, the abundance of one category usually prevailed, but the species richness was highest in different categories. This applies to number of individuals and species richness in the size categories B and C, species preferring open biotopes and species with low habitat preferences (eurytopic), and hygrophilous and eurytopic species (humidity affinity). The abundances of spring and autumn breeders were very similar, while the numbers of species differed substantially. A similar situation was found for the carnivorous and omnivorous species (Table 1).
The Simpson dominance and Berger-Parker indices were highest in locality SB3, where 80% of individuals were identified as Pterostichus melanarius. The second highest index was found in WS with a dominance of Pseudoophonus rufipes (70%), and the third highest in SB2, where Poecilus cupreus represented 60% of individuals. The dominance of P. rufipes in CB (47%) and of P. melanarius in SB1 (32%) was less pronounced. This was reflected in the species evenness, which was highest in SB1 and lowest in SB3 in which species evenness was very similar to WS. The Margalef index detected highest ratio between the number of species and the abundance in SB2 (59 species, 22,015 individuals), followed by SB1 (35 species, 5484 individuals), CB (35 species, 5831 individuals), WS (34 species, 9401 individuals), and SB3 (34 species, 15,573 individuals, Table 2).
The Jaccard and Sørensen-Dice indices showed dissimilarity between communities when all localities were compared. The most similar communities were found in the three geographically closest localities in South Bohemia (Table 3).

The Effect of Locality, Year, and Sampling Date on Carabid Communities (All Data Included)
The location explained 16.3% and 23.8% variability in the distribution of carabid species and functional categories, respectively. In the species-based canonical correspondence analysis (CCA), the WS locality The year explained 1.0% (F = 14.9, p = 0.001) and 1.9% (F = 27.8, p = 0.001) of variability in the species-and trait-based CCA, respectively. In the species-based CCA, each of the time series S and A explained 3.6% of variability (S: F = 56.3, p = 0.001, A: F = 54.7, p = 0.001). The joint analysis of both series proved their close correlation. In the trait-based CCA, the time series S explained 4.3% (F = 67.1, p = 0.001), and the time series A explained 4.7% of variability (F = 70.4, p = 0.001). The joint analysis of both series showed they were correlated, sharing 4.1% of the variability they explained. Only time series A was used in subsequent modeling.

Three Possible Ways of Using Carabids in PMEM
Variability explained by localities remained relatively stable in the SaS and SaT models (see Section 5.3) until species with abundance lower than 150 and 600 individuals, respectively, were disregarded ( Figure 1a). Subsequent step-wise elimination of the more abundant species resulted in a steep increase in explained variability (the difference between variability explained by two adjacent points in the graphs significantly increased, SaS: F 1,25 = 11.69, p = 0.002, SaT: F 1,25 = 35.84, p < 10 −5 ). Only a small increase in explained variability was observed with the TaT model ( Figure 1a). The curve derived from the SaS model intersects the curves of the TaT and SaT models at the points corresponding to species represented by 117 and 229 individuals, respectively, where the variability explained by the SaS model exceeded the variability explained by both the SaT and TaT models.  Variability explained by localities was about four-times higher than that explained by the sampling date (time series A) in the analysis that included all species. However, variability explained by the sampling date increased faster with species elimination, and eventually became half that explained by the localities (Figure 1b). The increase in variability explained by time series A was similar in the SaS and SaT models. The smallest increase in variability was observed in the TaT model ( Figure 1b). Variability explained by years was very low in all models (Figure 1c).
The lowest percentage of variability explained by localities occurred in the TaT model and in the SaS model at the beginning of modeling (x ≤ 150), suggesting that these two approaches were most appropriate to compare the least locality-dependent environmental impacts. However, these procedures required the determination of all individuals to the species level, and in the case of the TaT model also their classification into categories of functional traits. Neither of these requirements can be fulfilled in routine practice. However, variability explained by localities in the SaT model was at x = 229, exceeded by variability explained by the SaS model, and then increased slowly up to x = 600 (Table 4). Species with abundance lower than 230 (in total for all localities) could be excluded from the analysis, while species with abundance higher than 600 (10 species, 14 categories) had to be included (this range is highlighted in Figure 1) to avoid a high increase of variability. Six species were found in all localities and their total abundance surpassed 600 ( Table 5). The SaT model showed conditions that have to be fulfilled for reliable comparison of different localities in routine practice. In the SaT model at a cut-off level of 600, most of the categories of functional traits were distributed around the center of the ordination diagram; their presence was similar in all localities. However, the incidence of hygrophilous species tended to be higher in WS, and species in body size category C were most common in SB1 ( Figure 2). Similar carabid groupings were found in SB1 and SB2, which were about 2-km apart, including 1 km of a forest. A similar species composition was also found in SB3 (Table 3) When we compared plots with GM events and plots treated with insecticides at a cut-off level of 600, the variability explained by these localities was 20.7%. It was lower than baseline (Figure 1a: SaT model, x (no. of individuals) = 600, y (explained variability) = 28.1%), indicating low probability of an impact of GM maize on the carabid groupings in these localities (Figure 3a). When GM events were compared with near-isogenic cultivars, variability among these plots in different localities was 28.9% (Figure 3b), which is still around the level of variability explained by different localities. A higher difference would indicate that the GM crop had an impact on the agroecosystem.  In the SaT model at a cut-off level of 600, most of the categories of functional traits were distributed around the center of the ordination diagram; their presence was similar in all localities. However, the incidence of hygrophilous species tended to be higher in WS, and species in body size category C were most common in SB1 ( Figure 2). Similar carabid groupings were found in SB1 and SB2, which were about 2-km apart, including 1 km of a forest. A similar species composition was also found in SB3 (Table 3), but the species abundance was different ( Figure 2).
When we compared plots with GM events and plots treated with insecticides at a cut-off level of 600, the variability explained by these localities was 20.7%. It was lower than baseline (Figure 1a: SaT model, x (no. of individuals) = 600, y (explained variability) = 28.1%), indicating low probability of an impact of GM maize on the carabid groupings in these localities (Figure 3a). When GM events were compared with near-isogenic cultivars, variability among these plots in different localities was 28.9% (Figure 3b), which is still around the level of variability explained by different localities. A higher difference would indicate that the GM crop had an impact on the agroecosystem.

Discussion
Development of a PMEM method that is applicable in a large geographic area at reasonable cost is very challenging. Such a method should be based on indicators that are relatively easy to monitor, occur in vast territories from spring to autumn, and are exposed to the products of transgenes in GM crops. Carabids fulfill these requirements, but the feasibility of their monitoring has not been sufficiently analyzed. Therefore, our study focused on the nature and size of data needed for reliable distinctions between carabid communities in different localities, years, and sampling dates. Carabid assemblages were analyzed in respect to the species composition and functional traits. This approach was preferred over the analysis of total abundance of the carabid family.

Carabid Communities in Maize Fields
To facilitate comparison with other studies, we used several indices to characterize carabid communities in different localities. The values of diversity indices were more or less in the range of those reported across Europe [17,32,36,37]. Similarity indices based on the qualitative species comparison declined with the distance between compared localities. Conversely, functional diversity was similar in remote localities; all important categories of functional traits were present in all localities (TaT model, Figure 1a).
The Berger-Parker index is an effective, simple tool for monitoring impaired biodiversity in soil ecosystems due to human disturbances [38]. Index values increase from undisturbed to disturbed areas. In the present study, values ranging from 0.32 (SB1) to 0.70 (WS) were typical for sites with agricultural management and indicate the prevalence of one species. Use of this index can facilitate interpretation of soil biodiversity patterns in the context of ecosystem management and conservation. The index of species evenness was low. It ranged from 0.38 (SB) to 0.59 (SB1), and reflected dominance by one species and very low abundance of others [39].
The dominance of a few species varied in their dependence on locality, and annual changes in environmental and anthropogenic factors [36,37,39]. The 16 most abundant species in our study are common in the fields of Central Europe [17,34,36,37,39,40], and some are also found in southern Europe [9,11,41], the United Kingdom [32], and Balkan [42]. This suggests similarity in agrocoenoses across a large area; usually 10-20 frequently occurring species rotate in the position of the most abundant species. However, only five of these species (C. fuscipes, H. affinis, P. cupreus, P. melanarius, and T. quadristriatus) were among the 10 most common species in all localities. According to literature sources [9,11,17,32,34,36,37,[39][40][41][42] and our experience, we suggest that it is appropriate to compare the abundance of a group of 10-20 common species based on their functional classification. A similar conclusion has been reached by other authors [11].

Variability Explained by Locality and Environmental Variables
Site location is by far the most important source of variability [8]. Locality is therefore important for making a baseline between background variability caused by other factors than the studied treatment. Several studies have compared trials conducted under different management regimes and used different statistical methods to define the sample size sufficient to detect impacts of GM crops on the abundance of arthropods [8,9,11,[43][44][45]. We used CCA with three types of gradual species reduction or categories of functional traits to identify the minimum sample size that would sufficiently minimize the locality-dependent variability. Modeling was based on species abundance, and representation of categories of functional traits (explained in Section 5.3 Data Analysis). When all data were included in the analysis (first point in Figure 1a), the percentage of variability explained by localities was lowest in the SaS model, which showed the importance of species of low abundance for the similarity assessments of distant localities. The SaT model showed that variability explained by localities remains relatively low when species with abundance lower than 230 individuals are neglected. Species can be further eliminated from the analysis up to a threshold of 600 individuals per species. Exclusion of species with higher abundance causes a very steep increase in variability among localities. We do not recommend analyses based on less than 10 most common species. A trait-based analysis with the most common species is a compromise that can be utilized for large scale PMEM (limitations of PMEM are discussed in [10]). Species with an abundance lower than 230 can be excluded from the analysis, while those with an abundance higher than 600 must be preserved to avoid a high increase of unwanted variability.
We observed baseline background variability for carabid grouping comparisons in five different localities. If the percentages of variability explained by the GM and non-GM treatments were higher than the baseline for these localities, an impact of the GM crop on the agroecosystem may be indicated. The baseline for similar localities is lower than for the less similar localities. Thus, it is necessary to determine the approximate baseline variability for the examined localities.
Maize phenology was not crucial for the carabid grouping (shared variability between time series S and A). We propose to follow a time series based on the calendar date rather than on the date of maize sowing, as this varies between localities and years. The day 1 April seems to be a reliable landmark for Central Europe and probably for most of Europe.
Differences of 14 days between corresponding samplings in different localities were not uncommon in our study. Given the low percentage of variability explained by the sampling date, the precise timing of samplings in different localities is not needed. Our conclusions are mostly based on four samplings per season. Since the abundance and species composition fluctuate during the season [36], a minimum of four samplings per season are required. A reduced sampling number may lead to a considerable loss in the capacity to detect differences [45]. Most studies have confirmed the highest carabid abundance around the time of maize flowering [11,46], although this is not the rule [10]. It is advisable to cover the first part of the season (until grain development, [21]).
Unlike the findings of previous studies [8,11], the effect of different sampling years was relatively low and can be neglected. It seems that direct comparison with a current non-GM crop baseline should be used if available, but reference can also be made to historical baseline data [7].

The Applicability of Our Findings for GS in PMEM of GM Maize
Recommendations based on the findings of our study should be taken into account when designing statistical comparisons in PMEM. If any unintended effect is observed, the recommended data analysis will help to determine whether the adverse effect is associated with the use of a GM crop or whether it is a consequence of other environmental factors [47].
Many authors have highlighted the importance of field size and the availability of non-crop habitats adjacent to the field [36,39]. Those factors should be taken into account when differences among localities are to be interpreted. Multivariate analysis, as used in our study, is a multidimensional tool that considers the effects of many variables, and is appropriate for evaluating and subtracting the effects of these covariables.
Although the European Food Safety Authority (EFSA) claims that GS is not necessarily crop or event specific [4], we are convinced that GS methodology could be the same for different crops, but the direct comparison of carabid assemblage in various crops species, and the effects of agricultural practice are scientifically not defensible because the species and abundance composition is largely affected by crop type even in one locality. This was clearly shown for GM and near-isogenic cultivars of maize, beet, and oil seed rape [32].

Conclusions
The conservation of natural enemies, which as an important component of biological control, can be accepted as the endpoint of PMEM. We propose that these are monitored by analyzing captured carabid species with respect to their abundance and ecological functions. Reliable results are obtained with commonly occurring and abundant species representing important functional categories. Field location is the main factor limiting the detection of changes caused in carabid assemblages by other factors, including the introduction of GM crops. The location effect is preserved at a tolerable level by including only species that occur in relatively high numbers in all examined sites. We demonstrate that, in the case of our model, the inclusion of 10 species, each represented by ≥600 individuals (total count from all sites) is essential, while 70 species with less than 230 individuals each could be excluded without a significant increase in variability observed. The geographically closest localities with similar environmental properties should be preferably compared to reduce differences. At least four samplings during the season are recommended on similar, but not necessarily the same, dates. A time series based on the calendar date can be followed, and data from different years can be combined. Data from independent carabid analyses in maize fields (and possibly in the plantations of some other crops) can be included in future analyses and further increase the precision of PMEM. The proposed method is a compromise that enables the detection of small but meaningful differences at maximally reduced labor costs.

Experimental Localities
Field trials were performed in localities designated as South Bohemia 1, 2, and 3 (SB1, SB2, and SB3); Central Bohemia (CB); and western Slovakia (WS). Basic features of all localities are summarized in Table 6; details on WS site are provided here. Carabid data obtained from the sites SB1, SB2, and CB have been published (see references in Table 6) and are not evaluated from the perspective of the present study.
The WS trial was performed in a field previously planted with winter wheat. In the first trial, three treatments were applied in four replicates (12 plots in total, 30 × 30 m each). In the second trial, two treatments were tested in 10 replicates (20 plots in total, 10 × 10 m each).

Data Analysis
We used the following ecological indices to compare the diversity of carabid communities in different localities: Berger-Parker index, Margalef index, Simpson dominance index, and Species evenness. The Jaccard index and Sørensen-Dice index were applied to assess the similarity of communities in different localities [55].
Carabid distribution was analyzed using multivariate analysis (Canoco software for Windows 4.5, Plant Research International, [56]). The analysis concerned the abundance of species and their placing in the functional trait categories for body size, habitat and humidity affinities, breeding period, and food specialization. The detected gradient length (4.9) in the detrended correspondence analysis (DCA: 0.001 attributed to each value, detrending by segments, log transformation: x' = log (x + 1), downweighting of rare species) of distribution trends and the characterization of data enabled us to use canonical correspondence analysis (CCA: 0.001 attributed to each value, log transformation, Hill's scaling). The effects of geographic localization of each locality (dummy variables), year (dummy variables), and sampling date were tested (environmental variables). A two-time series was used to test the effect of sampling date: (1) time series, S: number of days from the sowing day, marked as number 1; and (2) time series A: number of days from 1 April, which was classed as Day 1. The day 1 April was selected based on the agro-technical term of maize sowing. The earliest possible term for sowing in Central Europe is around 5 April [57]. The joint explanatory effect of these two variables was assessed by the analysis of variability explained with the time series S and A in two separate CCAs, and together in a single CCA (variance partitioning procedure, [56]). The significance of the effects of environmental variables was tested by subtracting the effect of covariables (CCA in partial shape, Monte Carlo permutation tests, MCPT: 999 permutations, unrestricted permutations, forward selection). Covariables are environmental variables whose influence is subtracted before that of variables of interest is investigated [56]. Covariables were those environmental variables mentioned above whose effect were not tested in certain CCAs.
Variability explained by environmental variables in CCA was compared for three different types of gradual elimination of individuals with a cut-off level of 10: 1.
SaS model: the least abundant species were eliminated and a CCA was performed at the species level (Table S4); 2.
SaT model: the least abundant species were eliminated and a CCA was performed at the level of functional traits (Table S5); and 3.
TaT model: the least frequent categories of functional traits were eliminated and a CCA was performed at the level of functional traits (Table S6).
In the SaS and SaT models, species elimination proceeded until three most common species remained. In the TaT model, the categories of functional traits were eliminated until the three most common categories remained. Variability explained by models was compared using the curve (two-phase exponential association, coefficient of correlation r, Graph Pad Prism 4.5, [58]). One-way ANOVA (F-tests accompanied by degrees of freedom and degrees of freedom of the error) was applied to compare differences in variability explained during modeling [59].
We defined baseline as a background variability that is caused by other factors (covariables in multivariate analysis) than the variables (GM vs. non-GM) whose effect is important for the purpose of the study. We hypothesize that when we know baseline, it is possible to distinguish between the background variability and variability caused by growing GM maize. The example of separation of effect of baseline from effect of GM maize is presented in Figure 3a,b. The carabid grouping in GM maize is compared there with carabid groupings in plots treated with insecticides (as they are applied in most of maize cropping systems) and plots with near-isogenic cultivar, respectively.
Supplementary Materials: The following are available online at www.mdpi.com/2072-6651/9/4/121/s1, Table S1: The sample dates of deployment of pitfall traps in the locality western Slovakia in 2014 and 2015,  Table S4: Data for analysis of explained variability in SaS model, Table S5: Data for analysis of explained variability in SaT model, Table S6: Data for analysis of explained variability in TaT model.