Field trials to evaluate the effects of transgenic cry1Ie maize on the community characteristics of arthropod natural enemies

Possible non-target effect of transgenic cry1Ie maize exerts on natural enemy community biodiversity in the field is unresolved. In the present study, a 2-yr comparison of transgenic cry1Ie maize (Event IE09S034, Bt maize) and its near isoline (Zong 31, non-Bt maize) on natural enemy community biodiversity were compared with whole plant inspections, pitfall traps and suction sampler. Natural enemy diversity indices (Shannon-Wiener’, Simpson’s and Pielou’s index) and abundance suggested there were no significant differences between the two types of maize. The only exceptions were the Pielou’s index for whole plant inspections in 2013 and abundance for pitfall traps in 2012, which were significantly higher in Bt maize than those of non-Bt maize. The main species of natural enemies were identical in Bt and non-Bt maize plots for each method and the three methods combined. For whole plant inspections, Bt maize had no time-dependent effect on the entire arthropod natural enemy community, and also no effect on community dissimilarities between Bt and non-Bt maize plots. These results suggested that despite the presence of a relatively minor difference in natural enemy communities between Bt and non-Bt maize, transgenic cry1Ie maize had little, if any, effect on natural enemy community biodiversity.

A number of field trials have been conducted to assess the potential negative impacts of Cry1Ab, Cry1Ac and Cry3Bb1 maize on natural enemies. Results suggested Cry1Ab maize had no significant effect on population dynamics of ladybird beetles (Coleoptera: Coccinellidae) 21 , and also was compatible with other natural enemies from the families Anthocoridae, Coccinellidae, and Araneae 22 . In addition, a 2-yr field study suggested Bt maize varieties expressing insecticidal proteins of Cry34Ab1, Cry35Ab1, Cry1F also had no adverse effects on arthropod food web properties 23 . Laboratory studies feeding Trichogramma ostriniae (Chen & Pang) (Hymenoptera: Trichogrammatidae) with Cry1Ab maize pollen suggested Cry1Ab maize had no adverse effects on longevity and fecundity of T. ostriniae 24 . Likewise, Cry3Bb1 maize had no acute or chronic fitness effects on Coleomegilla maculata (DeGeer) (Coleoptera: Coccinellidae) that fed on aphid 25 . However, possible non-target effects on natural enemies also must be assessed for new types of Bt maize 26 . Such an assessment includes numerous factors, such as expression rates of the Cry protein in different plant tissues over the growing season, ingestability and susceptibility of natural enemies to the Cry proteins 12 , transfer probabilities of Cry proteins to higher trophic levels 17 , feeding ecology of both herbivores and natural enemies 27 , and influence of Cry proteins on predator behavior 28 .
The cry1Ie gene was first successfully identified from B. thuringiensis isolate Btc007 in Institute of Plant Protection, Chinese Academy of Agricultural Sciences 29 . Cry1Ie gene encodes insecticidal proteins toxic to ACB and cotton bollworm, Helicoverpa armigera (Hübner) (Lepidoptera: Noctuidae) 30 . This gene shows no cross resistance with Cry1Ab, Cry1Ac, Cry1Ah and Cry1F insecticidal proteins [31][32][33][34] . A study on environmental risk assessment of transgenic cry1Ie gene maize on arthropod biodiversity was conducted 35 , but did not include a focused analysis of the natural enemy's data. Because of the importance of natural enemies to the maize ecosystem, the present study will reevaluate the natural enemy data from that study in more detail.
The aim is to compare arthropod natural enemies in transgenic cry1Ie maize plots with those in near isoline Zong 31 maize plots with emphases on (i) species diversity and abundance, (ii) time-dependent effects of Bt maize on community composition, (iii) similarities of community structures between Bt and non-Bt maize and their response to maize type and sampling time. With this purpose, the abundance and diversity of natural enemies from the two-year field study were evaluated. The three multivariate techniques, redundancy analysis (RDA), principal response curve (PRC), and nonmetric multidimensional scaling (nMDS) were used to assess the community data of natural enemies in Bt and non-Bt maize plots.
We compared these significant diversity indices for sampling times with a two tailed t-test. All the diversity indices showed similar trends of temporal dynamics and had no consistent difference between Bt and non-Bt maize during the two years ( Fig. 2a-i). For each sampling time, the results of pairwise comparison showed that for whole plant inspections, Shannon-Wiener diversity index of Bt maize was significantly different from non-Bt maize at V12 stage in 2012 (t-value = 3.33, P = 0.029), and at V6 stage (t-value = 3.43, P = 0.024) in 2013 (Fig. 2a); Simpson's diversity index was significantly different at V12 stage in 2012 (t-value = 3.21, P = 0.033), and at V6 stage (t-value = 3.44, P = 0.026), R5 stage (t-value = − 2.79, P = 0.050) and R6 stage (t-value = − 4.04, P = 0.016) in 2013 (Fig. 2d); Pielou's evenness index of Bt maize was significantly greater than non-Bt maize at R6 stage (t-value = − 5.30, P = 0.013) (Fig. 2g), which may have accounted for the differences of Pielou's evenness index between Bt and non-Bt maize in 2013. No significant difference was found in the pairwise comparisons of the diversity indices between Bt and non-Bt maize in pitfall traps and suction sampler. Impacts of maize type and sampling time on the natural enemy abundance. During the two years, there were no significant differences in natural enemy abundance between Bt and non-Bt maize plots. The only exception was that in pitfall traps in 2012, natural enemy abundance in Bt maize was significantly higher The Y-axis shows the percentage of each taxa and the X-axis shows the maize type. This diagram illustrates dominant natural enemies (proportion > 1%) and the total proportion of taxa with percentages below 1% (others). Numbers above the columns show the total numbers of taxa collected with each method and three methods combined.
For each sampling method, natural enemy abundance showed similar trends of temporal dynamics ( Fig. 2j-l). No significant difference was found in the pairwise comparisons of natural enemy abundance between Bt and non-Bt maize at each sampling time.
Time-dependent effects of Bt maize on natural enemy community. Redundancy analysis (RDA) was performed to discern the possible relationship between natural enemy community composition and maize type, as well as sampling time. The results indicated that maize type and sampling time in total explained 12.72% of the variance of the natural enemy community data in 2012 (F = 4.15, P = 0.001; 999 Monte Carlo permutations), and 12.25% in 2013 (F = 3.98, P = 0.001; 999 Monte Carlo permutations). Among them, 11.91 and 0.81% of the variance in 2012, and 11.15 and 1.10% in 2013 were explained by the first and second axes, respectively (Table S2). When maize type and sampling time variations analyzed for contributions to natural enemy composition, maize type variance contribution (in 2012, R 2 = 0.11, P = 0.028; in 2013, R 2 = 0.14, P = 0.013; 999 Monte Carlo permutations) was significant, but appeared lower than that of the sampling time variance contribution (in 2012, R 2 = 0.78, P = 0.001; in 2013, R 2 = 0.68, P = 0.001; 999 Monte Carlo permutations) (Table S3).
Among all the canonical axes, only the first canonical axis was significant in 2012 (F = 7.78, P = 0.001; 999 Monte Carlo permutations) and 2013 (F = 7.25, P = 0.001; 999 Monte Carlo permutations), which explained 11.91 and 11.15% in 2012 and 2013 of the variation of the natural enemy community composition, respectively (Table S2). Principle response curve (PRC) analysis examined the time-dependent effects of Bt maize on the entire arthropod natural enemy community composition and considered the variations of non-Bt maize as baseline. PRC revealed that in 2012, 46.35% of the total variance of the natural enemy community data was explained by sampling time and 7.01% by maize type, as well as in 2013, 38.49% of the total variance was explained by sampling time and 10.20% by maize type. The statistical significance of the first PRC axes revealed no significant difference between the natural enemy communities in Bt and non-Bt maize plots in 2012 (F = 2.05, P = 0.969; 999 Monte Carlo permutations) (Fig. 3, Table S4), as well as in 2013 (F = 2.87, P = 0.597; 999 Monte Carlo permutations) (Fig. 3, Table S4). For simplicity Fig. 3 only summarized the species/families with species weights higher than 0.5 and lower than − 0.5 in each year. In 2012, M. tricuspidatus, and Syrphidae in Bt maize plots were higher than those in non-Bt maize plots, whereas C. sinica, Paederus Fuscipes (Curtis) (Coleoptera: Staphylinidae), Macrocentrus cingulum (Brischke) (Hymenoptera: Braconidae) and Clubionidae were more likely to follow the opposite trend. Among them, M. tricuspidatus (b k = 1.0124), and C. sinica (b k = − 0.9389) were likely to contribute most to the community response. In 2013, the abundance of Syrphidae, C. sinica, N. doenitzi and L. sinensis were higher in Bt maize plots than that in non-Bt maize plots, while the abundance of P. japonica, H. axyridis, Aphidiidae and M. cingulum were slightly lower in Bt maize plots (Table S5).

Diversity indices
Year Non-Bt maize Bt maize   Natural enemy assemblages in Bt and non-Bt maize plots. Nonmetric multidimensional scaling (nMDS) analysis (which clusters samples based on natural enemy composition) was used to visualize clusters of samples containing highly similar natural enemy composition (Fig. 4).
The nMDS plot showed that natural enemy communities in Bt and non-Bt maize were not well separated from each other for each sampling time during the two years (Fig. 4). This was supported by ANOSIM (in 2012, R 2 = 0.00, P = 1.000; in 2013, R 2 = 0.00, P = 0.908; 999 Monte Carlo permutations) (Fig. 4, Table S6). However  (Table 3).

Discussion
Most studies that have considered possible effects of Bt maize on non-target arthropods, including natural enemies, have focused on maize expressing Cry3Bb1 36 or Cry1Ac 37 insecticidal proteins. Only one study reported on possible effects of transgenic cry1Ie maize on the entire non-target arthropod biodiversity in the field 35 , and   no effects were found. Because of the ecological importance of natural enemies, the natural enemy subset of this study that included whole plant inspections, pitfall traps and suction sampler was evaluated in more detail.
Methods for analyzing the effects of environmental factors on species and ecological communities are generally fallen into two categories: (i) those that emphasize on the distributions of individual species, and (ii) those that emphasize on differences in the community composition 38 . Since one single method is not available to completely evaluate the natural enemy community, a combination of methods was necessarily performed to obtain a more comprehensive view of the natural enemy community in Bt and non-Bt maize plots. Therefore, diversity indices (diversity: Shannon-Weiner's, and Simpson's; evenness: Pielou's), abundance index, redundancy analyses (RDA), principal response curves (PRC) and nonmetric multidimensional scaling (nMDS) were used within this study.
The two diversity indices (Shannon-Weiner's, and Simpson's) and an evenness index (Pielou's) are useful indicators for measuring the disturbance of natural enemy communities 37 . Significant effect of maize type was observed only for Pielou's evenness index with whole plant inspections in 2013, where Bt maize was higher than that of non-Bt maize (Table 1). These results overall suggested that Cry1Ie Bt maize had little to no impact on natural enemy community biodiversity. However, when these indices were used to evaluate different stages of maize development, some differences were detected for whole plant inspections, but these differences showed no clear pattern (Fig. 2). Significant differences between Bt and non-Bt maize plots were observed at V12 stage in 2012 and V6 stage in 2013 for Shannon-winner index, at V12 stage in 2012 and V6, R5 and R6 stages in 2013 for Simpson's diversity index, and R6 stage in 2013 for Pielou's evenness index by whole plant inspections (Fig. 2a,d,g). Interestingly, in 2013 indices were higher in Bt maize plots than non-Bt maize plots for Simpson's diversity index during R5 and R6 stages, and Pielou's evenness index during R6 stage (Fig. 2d,g). Overall, these results suggested that the diversity indices between Bt and non-Bt maize had relatively minor differences and these differences showed no consistent pattern. In some cases, the Shannon-Weiner's and Simpson's diversity indices may have been influenced by rarely occurring species, which occurred more often in non-Bt maize compared to Bt maize, e.g., M. cingulum. No significant differences existed in natural enemy abundance between Bt and non-Bt maize except for pitfall traps in 2012, where we found fewer natural enemies in non-Bt maize plots compared with Bt maize plots ( Table 2). These findings were consistent with a previous 2-yr monitoring survey report that Bt maize expressing Cry1Ac proteins did not affect the natural enemy community 37 . Similarly, a 3-yr study reported that no detrimental effect of farm-scale Bt maize was observed on abundance of predatory arthropods in Spain 21 .
PRC is a multivariate technique used to assess the structure of arthropod and soil nematode community in Bt and non-Bt crop 37,39 . Two fundamental questions for Bt maize environmental risk assessment for natural enemies were answered by using PRC models for whole plant inspections in 2012 and 2013. First, does Bt maize alter the natural enemy composition in a series of repeated observation? Results clearly suggested there were no consistent differences over time as there were no significant effects of Bt maize on natural enemy population distribution over time when compared to non-Bt maize in 2012 and 2013 (Fig. 3). Second, how do individual natural enemy species respond to Bt and non-Bt maize (vertical line of right side of line charts; Fig. 3)? Among the 14 species/ families monitored in 2012, C. sinica, P. fuscipes, M. cingulum and Clubionidae were more abundant in non-Bt maize plots, but M. tricuspidatus and Syrphidae were more abundant in Bt maize. Among 22 species/families in 2013, P. japonica, H.axyridis, Aphidiidae and M. cingulum were more abundant in non-Bt maize, but Syrphidae, C. sinica, N. doenitzi, and L. sinensis were more abundant in Bt maize. The only consistent results were that Syrphidae was higher in Bt maize and M. cingulum, a specialist larval parasitoid, was higher in non-Bt maize. The latter is probably attributed to fewer Asian corn borers, O. furnacalis, in Bt maize.
The RDA allowed one to assess the relative contributions of maize type, sampling time and unknown factors with the abundances and species of natural enemies. Maize type and sampling time explained 12.72 and 12.25% of the variance in the natural enemy communities in 2012 and 2013, respectively. More than 80% of the variance was caused by other factors. Many factors may explain the differences besides genetic modification: (1) plot size and isolation among them determine the abundance of aboveground arthropods 40 ; (2) fewer target insects in Bt crop plots cause less damage in the plots, so the effects of herbivore-induced volatile in plants to natural enemies may be weak 41,42 ; (3) effects of the reduced number of target pests on natural enemy community abundance, and (4) transgenic crops with new genes could alter the physiological parameters 43 (i.e., lower N accumulation or differential plant development rate 44 ). A decline in target hosts is a likely mechanism to explain the reduction of  Table 3. Effects of maize type and sampling time on beta diversity (Bray-Curtis distance) of natural enemy community. *** P < 0.001.
Scientific RepoRts | 6:22102 | DOI: 10.1038/srep22102 natural enemies 45 . For example, fewer maize borers probably contributed to fewer M. cingulum, an Ostrinia larval parasitoid, in Bt maize plots. The Bray-Curtis distances for 2012 and 2013 showed that sampling time was a much more important factor than maize type in explaining natural enemy variability (Fig. 4). NMDS is less sensitive than the diversity indices to rare species occurrence in samples so may be more appropriate for comparing community compositions of Bt and non-Bt maize 46 .
Natural enemies contribute to the control of pest populations 47 , thus natural enemy biodiversity is linked to biological control of insect pests 48 . Transgenic crops can reduce pesticide applications 18 . If the abundance and diversity of natural enemies is maintained in transgenic crops because fewer insecticides are used, secondary pests are less likely to outbreak. Ultimately, Bt maize could foster durable integrated pest management 26,49 .
Assessment of possible impacts of Bt maize on natural enemy community and biodiversity is necessary for a pre-release environmental risk assessment 50 and post-release monitoring 51 . Taken as a whole, our study confirmed the results of similar field trials of Bt maize expressing Cry1Ac insecticidal proteins 37 demonstrating that Cry proteins expression in maize did not affect natural enemies populations (predators, parasitoids). Laboratory studies further confirm that direct feeding on Bt plant material poses a negligible risk for predators 22,52 . Field conditions are complex, and many factors cannot be controlled, weather, different prey types, food shortage, prey choice and competition between natural enemies and all these factors may influence the final results.
This study confirmed that despite relatively minor differences in natural enemy communities between Bt and non-Bt maize, transgenic cry1Ie maize had little effect on natural enemy community biodiversity. This 2-yr field assessment provided an assessment of possible ecological effects of Bt maize on natural enemy biodiversity. The indirect effects of Bt maize on natural enemy, such as the attraction of herbivore-induced volatile emission in maize on natural enemy 41 , or the abundance of food resources for natural enemies, were not separated from the direct effects. In addition, which natural enemies could ingest Bt insecticidal proteins when living/feeding on the hosts/prey were not determined. Such factors could be deciphered with additional long-term studies.

Methods
Ethics statement. For natural enemies sampling in maize field, no specific permits were required.
None of the species used in this study are endangered or protected. Sample collection. Whole plant inspections were carried out periodically for all the main growing stages of maize: 3 rd (V3), 6 th (V6), 9 th (V9) 12 th (V12) leaf stages, tasseling (VT), silking (R1), blister (R2), milk stage (R3), dent (R5) and physiological maturity (R6) in 2012 and 2013. In each Bt maize and non-Bt maize plots, twenty plants were randomly selected along two corner-to-corner diagonals (X shaped) (100 total plants per plot sampled each year). The numbers of visible natural enemies on stalks, leave, sheaths, tassels, husks and ears were quickly counted. If a species could not be identified, it was kept within a 5-ml plastic tube with 75% alcohol for later identification.
Pitfall traps and suction sampler were also used to collect natural enemies at V6, V9, R1 and R2 stages of maize development. Five sites were chosen in each plot, one site was in the plot center and the other four were in the middle of lines that connected the plot center with plot corners. Three pitfall traps spaced 0.5m apart in a line were established at each of the five locations in a plot. For each trap a plastic outer cup (15 cm in diameter × 10 cm depth) was buried in the ground with the upper rim of the cup level with the ground. Then one collection cup (5 cm in bottom diameter × 8 cm depth) with 75% ethanol, sugar, vinegar and water was put into each outer cup. Traps were exposed for about 24 h, after which they were collected and trapped natural enemies were stored in 75% ethanol until they were identified to species and counted. A Univac suction sampler (Burkard Manufacturing Co. Ltd, UK) was used to collect natural enemies, especially for the small individuals not easy to observe from 10 plants for each of the trapping locations. All the collected arthropods again, if possible, were identified to species.

Statistical analysis.
To analyze the natural enemy abundance and diversity, the number of species (species richness) and the relative abundances of individuals within each species (species abundance) per plot at each sampling time were calculated. The total number of individuals caught per plot at each sampling time was used as an index of relative abundance 53  Pielou's evenness index (J) was calculated as follows: where S is the total number of genera, H′ is Shannon-Weaver diversity index 55 .
Simpson's diversity index (D) was calculated as follows: N i = the number of individuals in the i th taxon; and N = the total number of individuals 56 . Those three diversity indices in whole plant inspections, pitfall traps and suction sampler in 2012 and 2013 were calculated using the 'diversity' function of the vegan package 57 in R. Total number of natural enemies were log(x + 1) transformed before calculating diversity indices.
Unequally spaced repeated-measures ANOVA ('proc mixed' procedure in SAS) 58 was used to test the effects of maize type and sampling time on natural enemy abundance, Shannon-Weiner's, Simpson's and Pielou's index. Maize type and sampling time were fixed factors, and block was a random factor. When a significant effect of sampling time was observed, comparison of mean values of abundance, H′ , J and D on each sampling time was conducted with a two-tailed t-test to detect significant differences between Bt and non-Bt maize plots.
One canonical analysis method (redundancy analysis, RDA), and its special type (principal response curve, PRC), and one indirect ordination method (nonmetric multidimensional scaling, nMDS) were used to compare differences in the composition of natural enemy communities and identify environmental controls of community composition statistically 38 . Here we only analyze the data collected by whole plant inspections in 2012 and 2013 due to the higher abundance of natural enemies. RDA was performed to investigate any linear relationships between the factors (maize type and sampling time) and the variation in natural enemy community composition. DCA was initially used to explore the range of variation within the natural enemy community composition 59 . If the longest gradient length was less than 3, the species were responding linearly to the gradients of explanatory variables, and therefore, linear response models were appropriate 60 . In our present study, the longest gradient length of DCA in 2012 and 2013 were short (1.74 and 1.52 SD, respectively), so RDA was chosen.
RDA is a constrained ordination technique in which the main axes are constrained to be linear combinations of the environmental variables 60 . Here we used RDA directly displayed the variation of natural enemy species as much as they can be explained by maize type and sampling time 61 using the 'rda' function of the vegan package 57 in R. The proportion of variance explained by sampling time and maize type was quantified using R 2 (adjusted) 62 . The total abundances of natural enemy were log(x + 1) transformed prior to analysis to stabilized variance and reduce the influence of dominant taxa on the ordination 63 . The level of significance of the canonical axes of RDA axes was tested by Monte-Carlo permutation test in R and yielded only one significant axis in each year. Thus the species-explanatory variables correlation along the first axis of RDA was used to set up the PRC.
PRC was used to visualize temporal changes in natural enemy composition caused by Bt maize as compared with non-Bt maize and also to quantify the contribution of each species to separate Bt maize from non-Bt maize (prc function, Vegan package in R). PRC is a special type of RDA for multivariate responses in a design with a series of repeated observation 64 . In the PRC analysis, principal component in species composition is plotted through time for explaining compositional deviations from the control represented as a horizontal line 65 . This is achieved by taking the non-Bt maize as the reference to Bt maize and by defining "time" as the horizontal axis of the diagram. The species weights (b k ) represent the highly affinity of the treatment, and indicate the direction of the changes in abundance 66 . Species weights between + 0.5 and − 0.5 were not displayed for their weak or unrelated responses to the PRC 67 . Monte Carlo permutations test was performed to test the significance between Bt maize and non-Bt maize for the axis of interest (in our case is the first axis), and the critical probability level for detecting significance between Bt and non-Bt was set at α = 0.01 67 .
The similarities of natural enemy community composition across maize type and sampling time were visualized using nMDS based on Bray-Curtis dissimilarity matrix 68 . A stress function ranged from 0 to 1 was used to assess the goodness of fit between the ordination and the original data of natural enemies (isoMDS function, Vegan package in R). The stress values were all below 0.2 (2012 was 0.1950, and 2013 was 0.1817), which suggested that the ordination was accurately represented the dissimilarity between samples 69 . Shepard diagram (stress plot function, Vegan Package in R) of non-metric fit illustrated that observed dissimilarities and the ordination distances were high correlated (non-metric fit was 0.962 in 2012, and 0.967 in 2013, respectively) (Fig. 6). In our study, the natural enemies surveyed in Bt and non-Bt maize by whole plant inspection in 2012 and 2013 distributed into 60 samples (10 sampling time × 2 maize types × 3 replications). Data were log (x + 1) transformed and normalized before Bray-Curtis dissimilarity were computed.
NMDS produced a two-dimensional graphical representation of the pairwise dissimilarity between natural enemy communities in each sample using the packages Vegan and Mass of R version 3.0.3 56 . Distances between each plot were estimated by using Bray-Curtis dissimilarity index 70 , which is one of the most robust statistics for multivariate ecological analysis of community composition among samples and slightly affected by the presence of rare species 46 .
NMDS also allowed us to determine which measured factors (maize type and sampling time) were significantly correlated with the nMDS ordination of natural enemy communities by returning squared correlation coefficients (envfit function, Vegan package in R) 57 . Fitted vectors were plotted onto the nMDS ordination of natural enemy community structures in Bt and non-Bt maize plots, and significance of fitted vectors were assessed by a permutation test 57 . We likewise used Bray-Curtis dissimilarity index to determine whether changes in natural enemy community composition between samples were related to changes in maize type and sampling time with ' Adonis' function in the vegan package 57 of R.
Statistical analyses were performed in the R software environment (v.3.0.3; R Development Core Team) 71 and SAS statistics package version 9.2 72 .