Impact assessment of genetically modified herbicide-tolerant cotton on arthropod communities

Cotton (Gossypium spp.) is one of the most important economic crops worldwide, and its production plays an important role in the economy of many countries. Genetically modified herbicide-tolerant (GMHT) crops, which were developed to minimize the losses caused by weeds, have gradually become the most widely adopted genetically modified crops in the world due to their economic and environmental benefits. However, the potential ecological and environmental risks of GMHT crops have attracted extensive attention and controversy. Arthropod communities form a prominent part of the biodiversity of agroecosystems and are important indicators of environmental health. Elucidating the effects of GMHT crops on the diversity of arthropod communities is necessary to ensure the safety of GMHT crops. In this 2-year study, we investigated the potential impact of GMHT crops on arthropod communities. The GMHT cotton variety GGK2 with glyphosate tolerance and its near-isogenic non-GMHT variety K312 were used for the experimental groups. The Shannon diversity index (H), Simpson diversity index (D), Pielou evenness index (J), and principal co-ordinates analysis (PCoA) of the Bray–Curtis distance were used to evaluate the population dynamics and biodiversity of arthropods in cotton fields. No significant differences were found between GGK2 and K312 in their total abundance of arthropod communities, and biodiversity indexes on most sampling dates. The arthropod composition in the GGK2 and K312 plots was similar. Sampling dates had a significant effect on biodiversity indexes, whereas no clear tendencies related to cotton variety or cotton variety × sampling dates interaction were recorded. In addition, PCoA revealed high similarity between the arthropod communities in the plots of the GMHT cotton variety GGK2 and its near-isogenic variety K312. There was no obvious difference in abundance, diversity indexes of arthropod communities between GMHT cotton variety GGK2 and its near-isogenic variety K312 under the small-scale planting regime.

Genetically modified herbicide-tolerant (GMHT) crops are a scientific breakthrough that have helped revolutionize weed management, greatly simplifying field weed control, reducing weeding costs, and increasing economic benefits (Duke 2015;García-Ruiz et al. 2018;Brookes and Barfoot 2020).
GMHT crops are the most widely adopted genetically modified crops in the world due to their economic and environmental benefits (Green 2012;Heap and Duke 2018). The planting area of GMHT crops has increased substantially since 1996, and the global planting area of GMHT crops reached 166.6 million hectares in 2019, accounting for 88% of the total planting area of biotech crops (ISAAA 2019). Although GMHT crops provide efficient tools for weed management, their potential ecological and environmental risks, such as gene flow, herbicide-resistant weeds, and potential biodiversity impacts, have attracted extensive attention and controversy all over the world (Li et al. 2012).
Arthropods play a prominent part in the biodiversity of the agroecosystem due to their species richness and many ecosystem functions, such as biological pest control, pollination, decomposition, and nutrient cycling. Moreover, they are important indicators of environmental health (Romeis et al. 2013;Schütte et al. 2017;Wang and Guan 2020). Thus, studying the effect of GMHT crops on arthropod communities is necessary (Li et al. 2012).
The risks of GMHT crops for arthropod communities have been reviewed in many publications, and most studies indicate that GMHT crops have no toxicity to arthropods and will not directly affect the diversity of arthropod communities (Buckelew et al. 2000;Bitzer et al. 2002;McPherson et al. 2003;Bohan et al. 2005;Jiang and Xiao 2010;Li et al. 2012).
In the current paper, we undertook a 2-year study with GMHT cotton variety GGK2 to further understand the effect of GMHT crops on arthropod communities. GGK2, which possesses the GR79 EPSPS/GAT gene, is a glyphosate-tolerant cotton variety developed in recent years. It shows five times more resistance to glyphosate with a ten-fold reduction in glyphosate residues. Thus, GGK2 has great potential for commercial adoption (Liang et al. 2017). However, only a few studies on the effects of cotton with the GR79 EPSPS/GAT gene on arthropod communities have been conducted.
Based on current knowledge, we hypothesized that compared with non-GMHT cotton, GGK2 does not affect arthropod communities in terms of their composition, abundance, and biodiversity indexes. In this study, the Shannon diversity index (H), the Simpson diversity index (D), the Pielou evenness index (J), and principal co-ordinates analysis (PCoA) based on the Bray-Curtis index were used to evaluate the population dynamics and biodiversity of arthropods in the field. This approach can help us clarify the potential effects of the GMHT cotton variety GGK2 on arthropod communities.

Cotton varieties
The GMHT cotton variety GGK2 with the GR79 EPSPS/ GAT gene and its near-isogenic non-GMHT variety K312 were used in this study. The seeds of these varieties were provided by Biotechnology Research Institute, Chinese Academy of Agricultural Sciences.

Experimental design
In 2019 and 2020, two cotton varieties, GGK2 and K312, were grown in an experimental field (36°05′03″N, 114°30′55″E) located in Anyang, Henan, China. The growing season for cotton is from April to October. Cotton was planted on May 3 and April 24 in 2019 and 2020, respectively. The distance between cotton plants was 25 cm, and the row spacing was 60 cm. A total of 5-10 cotton seeds were sown for one plant seedling. After sowing, cultivation management was carried out in accordance with the local conventional cultivation management mode. Normal irrigation, thinning, seedling replenishment, pruning, topping, and weeding practices were performed. The cotton plants were not treated with any pesticides during the field experiment. The experiment was carried out in a randomized complete block design with two varieties replicated three times. Each plot had dimensions of 10 m × 15 m.

Sampling design
Arthropod numbers were observed and recorded throughout the whole cotton growing season through direct visual examination with some modifications. The sampling time was from 8:00 AM to 12:00 AM. In each plot, 25 cotton plants were sampled by following an "X" pattern that covered the whole plot. The investigation was performed by sampling five cotton plants in each location, which contained four corners and a center. The sampling was started at approximately 2 m into the plot to avoid edge effects. The number and species of all arthropods above ground on each plant were determined by careful searching on the cotton stems and both sides of the leaves. Active insects were identified first. Unknown species were collected and preserved in 75% alcohol for later identification in laboratory. Sampling took place every 7 days from seedling emergence to harvest. The investigation was conducted consecutively from mid-May to late September.
In accordance with the method of Wang and Guan (2020), the arthropod communities were roughly divided into four groups: predatory, herbivorous, parasitic, and neutral (including saprophytic groups and groups that are neither harmful to plants nor eat other insects). Species in Araneida (Arachnida) were classified as predatory.

Statistical analysis
Three indexes were used to analyze the diversity of arthropod communities: the Shannon diversity index (H), the Pielou evenness index (J), and the Simpson diversity index (D). These indexes were calculated by using the following formulas, where P i representing the ratio of the number of the ith individual to the total number of individuals, S representing the total number of a species, N i representing the number belonging to the ith type and N representing the total number of arthropods: The abundance of arthropod communities, H, D, and J were analyzed by using repeated measures ANOVA with block as the random effect, varieties as the fixed effects, and sampling dates as the repeated-measures. The abundance of arthropod communities was subjected to log 10 (x + 1) transformation before analysis. The mean values of H, J, and D for each sampling date were compared by using Student's t test to detect significant differences between GGK2 and K312.
In this study, H, D, and J were calculated by using the diversity function of the "vegan" package in R (Oksanen et al. 2020). The P values of H, J, D, and the abundance of the main groups were calculated by using SPSS 25.0 (IBM Corp., Armonk, USA). Images were generated by R 4.0.3 (The R Core Team 2020) and GraphPad Prism 6.0 (GraphPad Software, San Diego, USA). The statistical threshold was 0.05 for all the tests.

Principal co-ordinates analysis (PCoA) of Bray-Curtis distance
The Bray-Curtis index is an important indicator to compare the differences between two samples in cluster analysis. The value ranges from 0 (similar) to 1 (different). Bray-Curtis distances between arthropod communities from each variety plot at each sampling date were calculated by using the vegdist function of the vegan package in R (Oksanen et al. 2020). PCoA based on Bray-Curtis distances was performed by using the pcoa function of the ape package in R (Paradis et al. 2019).

Arthropod communities in the GGK2 and K312
A total of 96 476 arthropod individuals (GGK2: 44 834, K312: 51 642) were obtained from the survey in 2019 and 84 412 arthropod individuals (GGK2: 42 325, K312: 42 087) in 2020. These individuals belong to 37 families of 13 orders. The most abundant families were Aleyrodidae, Aphididae, Thripidae, Psocoptera, Linyphiidae, Anthocoridae. In this 2-year study, significant differences between GGK2 and K312 were observed in the abundance of Araneae (predatory) (P = 0.017) and Coleoptera (predatory) (P = 0.037) in 2019 and Trombidiformes (herbivorous) (P = 0.041) and Coleoptera (herbivorous) (P = 0.046) in 2020. However, no statistically significant differences were observed between GGK2 and K312 in the individual numbers of arthropods and predatory, parasitic, herbivorous, and neutral groups (Table 1). We analyzed the arthropod composition of GGK2 and K312 at the order level to investigate the influence of GGK2 on arthropod composition. Our results showed that the arthropod composition in GGK2 and K312 plots was identical (Fig. 1), and the main orders were Hemiptera and Thysanoptera. These two orders accounted for more than 80% of the total abundance of arthropod communities (Fig. 1). Further statistical analysis was performed on the individual number of arthropods. No significant differences were observed between GGK2 and K312 in the abundance of the main orders over the 2-year study except for Araneae (P = 0.016) in 2019 and Trombidiformes in 2020 (P = 0.041) (Fig. 2).

Effects of GMHT cotton on the dynamics of arthropod communities in the field
An investigation was conducted consecutively from mid-May to late September to understand the trend of fluctuation in community structure. We performed 17 investigations each year. The fluctuation of arthropod communities' structure is shown in Fig. 3. The community composition of arthropods in the GGK2 and K312 plots changed over time throughout the entire 2019 and 2020 cotton growing season. Hemiptera and Thysanoptera were the two orders with the largest proportion throughout the whole cotton growing season. The proportion of Hemiptera insects decreased first and then increased, whereas that of Thysanoptera insects increased first and then decreased. The proportion of arthropods in other orders was relatively small, and also changed with time.

Effects of GMHT cotton on the biodiversity indexes (H, D, and J) of arthropod communities
We analyzed the biodiversity indexes of the arthropod communities to further understand the differences between the two cotton varieties. The results showed that GGK2 and K312 had no significant differences in H, D, J (P = 0.133, 0.050, and 0.339 for H, D, and J, respectively) in 2019; the same results were obtained in 2020 (P = 0.580, 0.655, and 0.614 for H, D, and J, respectively) ( Table 2 and Fig. 4). In addition, during   (Table 2). We analyzed the diversity indexes of arthropod communities throughout the whole growing season to clarify the variation of diversity indexes over time. In 2019, H, D, and J exhibited a similar rise-decline-stabilize trend. In 2020, H, D, and J all showed consistent fluctuations with peaks at weeks 2-3, 6-7, and 12, respectively.
The sampling dates with significant differences between GGK2 and K312 were the second week of 2019 for D (P = 0.039); the 15th week of 2019 for H (P = 0.008), D (P = 0.002), and J (P = 0.044); the first week of 2020 for H (P = 0.026) and D (P = 0.014); the second week of 2020 for D (P = 0.026); and the ninth week of 2020 for H (P = 0.032) (Fig. 5).
In summary, our results showed no statistically significant differences between GGK2 and K312 in H, D, or J on most sampling dates, although some differences between GGK2 and K312 were obtained on a few sampling dates.  In addition, H, D, and J showed different fluctuation trends between 2019 and 2020 (Fig. 5).

Principal co-ordinates analysis (PCoA) of Bray-Curtis distance of arthropod communities of GMHT cotton and control
The results of principal co-ordinate analysis (PCoA) based on Bray-Curtis distance between arthropod communities from each variety at each sampling date are shown in Fig. 6. The samples of GGK2 and K312 collected on the same sampling date were clustered together. However, the samples fluctuated dramatically over sampling dates. This pattern was consistent within the two varieties over 2 years. This consistency indicated that cotton growth time and developmental stage were the main factors influencing arthropod communities.

Discussion
GMHT crops are the most widely used biotech crops in the world and play an important role in simplifying farming measures and reducing weed management costs (Duke 2005;Gianessi 2008). However, their potential ecological and environmental risks have attracted extensive attention (Li et al. 2012). Arthropods are the most diverse animal species in the world and play a prominent part in agroecosystems. Thus, the investigation of arthropod communities has become an important step in environmental safety assessment of genetically modified crops (Romeis and Widmer 2020). Unlike genetically modified insect resistant Bt crops, GMHT crops do not produce insecticidal proteins and thus they should have no direct toxic effects on arthropods (Lu 2008). At present, most studies have shown that GMHT crops have no direct impact on arthropod communities. Jiang et al. reported that the GMHT rice Bar68-1 had no significantly affect on the composition and diversity of arthropod communities in the canopy (Jiang and Xiao 2010). Similarly, Liu et al. reported that there was no significant difference in arthropod diversity (the number of insects, H, D, and J) between GMHT soybean ZUTS-33 and non-genetically modified control soybean HC-3 and ZH-13 (Liu et al. 2020). In addition, the ecological risk assessment results of GMHT maize showed that the arthropod communities structures in the fields of the transgenic maize C0010.1.1, its counterpart 178, and the conventional line Xianyu 335 are equivalent, and found that the genetically modified maize had no significant effects on field arthropod communities (He et al. 2018).
In this work, the effects of GMHT crops on arthropod communities were evaluated with the GMHT cotton variety GGK2 and its near-isogenic non-GMHT variety K312. No significant difference was found between GGK2 and K312 in the total abundance, biodiversity indexes of arthropod communities on most sampling dates. Arthropod composition was identical in GGK2 and K312 plots. Sampling dates had a significant effect on biodiversity indexes, but no clear tendencies related to cotton variety or cotton variety × sampling dates interaction were recorded. In addition, PCoA results indicated the high similarity of arthropod communities in the plots of the GMHT cotton variety GGK2 and its near-isogenic variety k312. In summary, our results are consistent with the findings of previous research.
GR79 EPSPS is a novel class II EPSPS that was identified in recent years, and the GMHT cotton variety GGK2 with GR79 EPSPS/GAT gene is a new cotton variety with glyphosate tolerance (Liang et al. 2017). To the best of our knowledge, only a few studies have assessed the safety of cotton with GR79 EPSPS/GAT gene on arthropod communities (Liang et al. 2017). Further research should be undertaken to assess the potential effects of GGK2 with glyphosate spraying on arthropod communities.
In general, GMHT crops are non-toxic to farmland animals and may not have a direct effect on the farmland ecosystem (Lu 2008). However, the largescale commercial planting of GMHT crops may cause changes in the farmland environment by changing weed species and management measures, and may potentially exert an indirect impact on the farmland ecosystem (Lu 2008;Li et al. 2012). Therefore, future studies should not only focus on the short-term impacts of GMHT crops on agroecosystems, but also accumulate additional relevant data to observe the long-term effects on agroecosystems.

Conclusion
Through the analysis of the abundance, composition, Shannon index (H), Simpson index (D), Pielou index (J) and principal co-ordinates analysis (PCoA), no obvious difference was found between the arthropod communities in the plots of GMHT cotton variety GGK2 and its near-isogenic variety K312 under the small-scale planting regime. There were no significant differences between GGK2 and K312 in abundance and biodiversity indexes on most sampling dates. PCoA revealed that there was high similarity between the arthropod communities in the plots of GGK2 and K312.