Variation among conventional cultivars could be used as a criterion for environmental safety assessment of Bt rice on nontarget arthropods

The current difficulty facing risk evaluations of Bacillus thuringiensis (Bt) crops on nontarget arthropods (NTAs) is the lack of criteria for determining what represents unacceptable risk. In this study, we investigated the biological parameters in the laboratory and field population abundance of Nilaparvata lugens (Hemiptera: Delphacidae) on two Bt rice lines and the non-Bt parent, together with 14 other conventional rice cultivars. Significant difference were found in nymphal duration and fecundity of N. lugens fed on Bt rice KMD2, as well as field population density on 12 October, compared with non-Bt parent. However, compared with the variation among conventional rice cultivars, the variation of each parameter between Bt rice and the non-Bt parent was much smaller, which can be easily seen from low-high bar graphs and also the coefficient of variation value (C.V). The variation among conventional cultivars is proposed to be used as a criterion for the safety assessment of Bt rice on NTAs, particularly when statistically significant differences in several parameters are found between Bt rice and its non-Bt parent. Coefficient of variation is suggested as a promising parameter for ecological risk judgement of IRGM rice on NTAs.

The current difficulty facing risk evaluations of Bacillus thuringiensis (Bt) crops on nontarget arthropods (NTAs) is the lack of criteria for determining what represents unacceptable risk. In this study, we investigated the biological parameters in the laboratory and field population abundance of Nilaparvata lugens (Hemiptera: Delphacidae) on two Bt rice lines and the non-Bt parent, together with 14 other conventional rice cultivars. Significant difference were found in nymphal duration and fecundity of N. lugens fed on Bt rice KMD2, as well as field population density on 12 October, compared with non-Bt parent. However, compared with the variation among conventional rice cultivars, the variation of each parameter between Bt rice and the non-Bt parent was much smaller, which can be easily seen from low-high bar graphs and also the coefficient of variation value (C.V). The variation among conventional cultivars is proposed to be used as a criterion for the safety assessment of Bt rice on NTAs, particularly when statistically significant differences in several parameters are found between Bt rice and its non-Bt parent. Coefficient of variation is suggested as a promising parameter for ecological risk judgement of IRGM rice on NTAs.
To meet the demand for food in the face of relatively limited arable land, China has devoted great efforts into developing genetically modified (GM) crops, especially insect-resistant GM (IRGM) rice lines. Cry proteins isolated from Bacillus thuringiensis (Bt) are the most widely used insecticidal proteins in IRGM rice. Since the first Bt rice plant was developed in 1989, over a dozen Bt rice lines with high resistance to lepidopteran target pests have been developed [1][2][3][4][5][6] . Two Bt rice lines, Bt Shanyou 63 and Huahui-1, received biosafety certificates in Hubei province in 2009, but neither has yet been approved for agricultural production. The issues related to the commercialization of GM crops include ecological risk, food safety, biosafety regulation, adoption by farmers and public acceptance. A relatively well-developed regulatory system for risk assessment and management of GM plants has been developed in China 7 . It was predicted that farmers would value the prospect of increased yields and the reduced use of pesticides and would readily adopt the production of Bt rice, based on experiences with Bt cotton and virus-resistant papaya 8,9 . The main factor slowing the pace of commercialization of GM rice in China is low public acceptance, which arises out of fear for human health and the environment 10 .
Food safety assessments of GM crops have been conducted investigating both their intended and unintended effects. Intended effect assessments have focused on measuring the thermal stability, digestibility, toxicity and allergenicity of introduced proteins as well as their metabolites. Unintended changes were assessed through compositional comparisons between transgenic and non-transformed parent plants following the principle of

Results
Biological parameters of N. lugens on different rice cultivars in the laboratory. Nymphal development duration. The nymphal development duration of N. lugens fed on Bt rice lines was approximately 18.5 days, while on their non-Bt parent it was 17.3 days. When analysed independently, the nymphal duration of N. lugens was longer when fed on Bt rice versus on the non-Bt parent XS11, especially for insects fed on KMD2 (Fig. 1, F = 3.4135, df = 2,163, p = 0.0428). The coefficient of variation (C.V) among Bt rice lines and non-transgenic parent was 6.3%. The range of N. lugens nymphal duration among conventional japonica rice cultivars was 16.5 to 19.0 days, with a mean at 17.7 days and coefficient of variation at 7.1%; and the variation of nymphal duration on conventional indica rice cultivars was even larger (15.0 to 24.0 days, C.V, 18.5%; Table 1).
Survival rate. The survival rates of N. lugens nymphs fed on Bt rice lines and the non-Bt parent, together with 14 other rice cultivars, are shown in Fig. 1 and Table 1. The coefficient of variation in survival rates among N. lugens nymphs fed on Bt rice and non-Bt parent was very small (7.0%), compared with that of insects fed on conventional japonica rice cultivars (13.4%) and indica rice cultivars (32.8%). No statistically significant difference was detected in survival rate of N. lugens fed on Bt rice lines (both KMD1 and KMD2) from that of insects fed on the non-Bt parent XS11 (F = 2.333, df = 2,17, p = 0.1780). Only N. lugens fed on IR72 and IR42, which contain the N. lugens resistance genes bph2 and Bph3, respectively, had significantly lower survival rates than all of the other treatments. Bt rice and non-Bt parent was high (52.0%), but it was still lower than the C.V among different conventional rice types (60.3% for japonica rice and 107.1% for indica rice). The honeydew produced by N. lugens female adults fed on rice cultivars carrying resistance genes (IR26, IR72, IR42) was much lower than the others. However, statistical difference was only seen when females fed on the hybrid indica rice ZZY1 ( Fig. 1 and Table 1).
Fecundity. The fecundity of N. lugens on KMD2 was significantly lower than that on XS11 when we compared Bt rice with the non-Bt parent independently (F = 3.6750, df = 2,44, p = 0.0405). However, mean levels of fecundity of females fed on the conventional japonica rice cultivars and hybrid indica rice cultivars were similar to those of insects fed on Bt rice lines ( Fig. 1, Table 1). Meanwhile, the C.V of fecundity among Bt rice and non-Bt parent (21.3%) was smaller than that among different conventional rice types (31.9% for japonica rice and 64.8% for indica rice). And the fecundity of N. lugens on Bt rice KMD2 (177.3 eggs per female) fell within the 95% confidence interval of conventional japonica rice cultivars. Only the N. lugens females fed on hybrid rice cultivar LYPJ laid significantly more eggs than those fed on all other rice cultivars except TN1 and ZZY1. By contrast, significantly fewer eggs were laid by N. lugens fed on IR42 than on most of the other rice lines.
Population density of N. lugens on different rice cultivars under field conditions. There was no significant difference in annual mean populations of N. lugens in the field among most of the cultivars, as shown in Table 2. No significant difference between Bt rice lines and the non-Bt parent was found when analysed independently either (F = 3.7641, df = 2, 8, p = 0.1204). The C.V of annual mean population between Bt and non-Bt parent was 21.8%, which was much smaller than that among conventional japonica or indica rice cultivars (58.4% and 67.6% respectively). Except for two rice lines, IR72 and ZJ22, on which the N. lugens population was consistently low throughout the experimental season, the population density varied significantly among sampling dates along with the development stages of N. lugens and rice (Table 2). Therefore, we further analysed the data from each sampling date. On 5 September, the C.V of N. lugens population density among Bt rice and its non-Bt parent (27.7%) was similar to that among conventional japonica cultivars (31.2%); but it was much smaller than those among conventional or hybrid indica rice cultivars. By contrast, on 24 September and 12 October, the C.V of N. lugens population density among hybrid indica rice was lower than that among conventional indica or japonica rice cultivars. Nevertheless, they were all higher than the C.V between Bt rice and non-Bt parent, although a significantly lower N. lugens population was detected on Bt rice KMD2 on 12 October when Bt rice lines were compared with the non-Bt parent XS11 independently (Supplementary Figure

Principal component analysis (PCA).
We conducted PCA to identify the contribution rate of each parameter to the performance of N. lugens in the laboratory and in field (Fig. 2). The results showed the first three eigenvalues corresponded to ca 87.9% of the accumulated contribution. All 17 samples were represented The conventional rice cultivars were devided into three groups as conventional japonica rice, hybrid indica rice and conventional indica rice. In each low-high bar graph, the left border represents minimum value in the category, while the right border represents the maximum value; the line in the box represents the mean. Statistical difference was tested only between Bt lines and the non-Bt parent. *Indicates a significant difference according to one-factor ANOVA analysis and Tukey's multiple-range test (p < 0.05).
two-dimensionally using their PC1, PC2 and PC3 scores in two separate plots. PC1 explained 54.7% of the variation and showed a separation of rice cultivars IR42 and IR72. The distribution of rice cultivars on PC1 was mainly determined by laboratory biological parameters, especially nymphal duration, survival rate and fecundity. PC2 accounted for 19.1% of the total variation and showed a separation of rice cultivars including IR26. The distribution of rice cultivars on PC2 was mostly affected by field population density. PC3 accounted for 14.1% of the variation. The weight of honeydew and field population density positively affected the distribution of rice cultivars on PC3. The two Bt rice lines did not show marked separation on PC1, PC2 and PC3 from their non-Bt parent XS11.

Discussion
In the present study, the biological parameters and field abundance of N. lugens on Bt rice lines KMD1/KMD2 were compared with those on the non-Bt parent XS11 and 14 other conventional rice cultivars. Compared to variations in rice lines developed through modern biotechnology from their non-transgenic parent, variations in conventionally bred rice cultivars are even larger; however, they are usually accepted by the public without hesitation. We propose the use of variation in NTA biological parameters derived from a set of conventional cultivars with a history of safe production as a criterion for safety assessment of Bt rice on NTAs. When statistical differences were detected between Bt rice and the non-transgenic parent, the variation of the parameters in question were compared with the variation range of commercial rice lines which are considered to be normal for the crop. Then, the question of whether Bt rice is as safe as conventional rice can be answered. It is similar to that used for food and feed risk assessments 12 . The German Advisory Council on the Environment also defined harm as changes that go beyond the natural range of variability for a particular asset of value 51 .
The Bt rice line KMD2, which was reported to prolong the nymphal duration and affect the fecundity of N. lugens 45 , poses no harm when environmental safety is a protection goal. However, there have been reports of positive effects on nontarget pests or negative effects of Bt crops on nontarget arthropods. Mirid bug (Heteroptera: Miridae) population sizes increased in cotton and multiple crops were correlated with wide-scale adoption of Bt cotton in China 53 . Bt11 and Mon810 maize showed remarkable positive effects on the performance of the corn leaf aphid Rhopalosiphum maidis (Hemiptera: Aphididae) 54 . Increased survival rate of N. cincticeps on Bt rice T2A-1, and longer larval developmental time of the predator P. japonica have been reported when pollen of Bt rice T2A-1 and T1C-19 were used as food 36,41 . Parasitoid mummies are less abundant in Bt cotton CCRI 41 plots compared to conventional cotton plots 55 . For such situations, our proposal of using variation as a guideline would be helpful in safety judgment.
In the current study, low-high bar graphs demonstrate visually the range of each biological parameter of N. lugens on different rice types. Larger variations in range can be seen in conventional rice cultivars, especially in   From the two score plots, we can also see clearly that the variation between Bt rice and the non-Bt parent is small. As shown in Fig. 2, only the two resistant rice cultivars IR42 and IR72 are distinctly different from others. And only IR26 is separated from others due to high annual population density but low honeydew weight. Both low-high bar graphs and PCA score plots provide evidence that, compared with the non-transgenic parent, the resistant level of Bt rice to N. lugens was not altered by transgenic manipulation. We calculated the 95% confidence intervals and coefficient of variation (C.V) of each parameter among different rice types to quantify the varation range. C.V is a standardized measure of dispersion of a probability distribution or frequency distribution. It shows the extent of variability in relation to the mean of the population, and is widely used as an index of reliability or variability in medical and biological sciences 56 . The confidence interval is also estimated based on the probability, but it is susceptible to sample size and distribution. Compared with the 95% confidence interval, the coefficient of variation might be a more promising parameter for such studies. In our case, the C.V between Bt rice and non-Bt parent was closer to that of conventional japonica rice or hybrid indica rice cultivars. Furthermore, most of the parameters investigated on Bt rice fell within the 95% confidence interval of conventional japonica rice or hybrid indica rice. Meanwhile, the background variability is high for field population densities especially among conventional indica rice cultivars, which was due to the inclusion of two N. lugens resistant rice lines (IR42, IR72). Although the field study was relatively small, due to logistical and regulatory constraints, the findings of the field study were supportive of the findings in laboratory study. As for most rice cultivars, the field population density represents both their attractiveness and tolerance. Inconsistent developmental stage of rice, sampling date and agronomic performance, such as greater tillers number, stronger stems and taller plants might all affect the results of field abundance investigation. So, it would be better to begin field investigation of N. lugens in late August with a sampling interval of 7-10 days. Rice cultivars with similar agronomic characters to the Bt rice evaluated should be used and a resistant and a sensitive comparator should also be set. Traditionally cultivated crops with a history of safe use for consumers/domesticated animals have already been used as comparators in food and feed risk assessments according to the guideline of EFSA (2011) 57 . The current report presents the range of variation of different rice on NTAs. According to our observations, both laboratory and field experiment revealed a larger variation range in biological parameters and field abundance of N. lugens among conventional rice cultivars. Our results help confirm the notion that the natural variation range could be used as a criterion for environmental risk evaluation of GM crops on NTAs. When significant differences are found on NTAs between GM crops and comparators, the question of whether it is as safe as conventional rice could be answered by comparing the C.V from GM crops and their comparators with those among conventional rice cultivas, especially for situations where significant effects on nontarget arthropods of IRGM crops were found. Only those that fell outside the normal variation range should be suggested for further evaluation in terms of safety 12 . Further experiments are needed to develop well-established models for practical use in risk assessment of Bt rice on NTAs.

Methods
Experimental materials. Two Bt rice lines (KMD1 and KMD2) developed from two T 0 plants expressing Cry1Ab driven by the maize ubiquitin promoter, as well as the untransformed parental cultivar Xiushui11 and 14 other non-transgenic rice cultivars, were used for laboratory and field evaluations. Both Bt rice lines had high resistance to stem borers and leaf folder under laboratory and field conditions 58 . The 14 non-GM rice cultivars included three conventional japonica rice cultivars, two conventional early season indica rice cultivars, three semilate indica rice cultivars with N. lugens resistance genes and six hybrid indica rice cultivars (Table 3).

Insects.
A colony of N. lugens was collected from the paddy field at the experimental farm of Zhejiang University, Hangzhou, Zhejiang Province, China, in 2011 and reared on susceptible 'Taichung Native1' (TN1) rice seedlings in nylon mesh cages in a phytotron (22 ± 2 °C, 60-70% relative humidity, and a 14:10 h light: dark regime.  For reproduction analysis, the sex of each adult was determined on emergence. A newly emerged female and male from the same tested rice plant were mated and introduced onto a 60-d-old plant of the same type. Each plant was maintained individually in a plastic bottle (10 cm in diameter and 30 cm in height) covered with nylon mesh at the top for ventilation. Five to 15 biological replicates were performed per rice line. The plants were changed every 5 d until the adults died. The number of eggs laid per female was individually recorded with the aid of a stereomicroscope. Honeydew was also collected from individual female adults feeding on 60 d-old plants of the same type through the Parafilm sachets. Parafilm sachets were prepared as described by Heinrichset al. 59 and weighed. A previously starved 2-d-old female adult was transferred into each sachet, and then wrapped around the rice stem carefully. After feeding for 24 h, we removed the N. lugens and weighed the sachets containing the honeydew again. The difference between the two weights was the weight of honeydew 59 . Six to 12 biological replicates were prepared per rice line. All experiments were conducted in the same phytotron where the N. lugens were reared. Field Experimental Design. A total of 17 rice lines, including two Bt rice lines and their non-transgenic parent as well as 14 other conventional rice cultivars, were used for field studies at Changxing Agrotechnical Experiment Farm of Zhejiang University, Zhejiang, China in 2011. Rice seeds were sown on 1 July and transplanted on 30 July. The experiment was performed with a randomized block design with three blocks × 17 rice lines. Therefore, the experimental field was divided into 51 small plots. Each experimental plot was 2 × 1.5 m in size and separated on all sides by a 30-cm-wide walkway. Seedlings were hand transplanted at a rate of three seedlings per hill spaced 16 × 16 cm apart. The entire experimental field was surrounded by four border rows of TN1. Normal cultural practices for rice cultivation were followed during the entire experimental periods except that no insecticides were applied. The densities of N. lugens adults and nymphs were sampled by the beating tray method as described by Chen et al. 45 . On each sampling date, 5 hills were sampled at random along a diagonal line in each plot. Field population density was investigated every 20 days throughout the season beginning at the tillering stage.
Statistical analysis. Data including nymphal survival rate, developmental duration, honeydew weight and fecundity were analysed by one-way ANOVA in a completely randomized design using the Data Processing System (DPS) package Version 15.10 60 , followed by Tukey's multiple-range test. Population densities of N. lugens in the field were analysed using the GLM model repeated-measures analysis of variance in SAS v.9.1, where date was used as repeated factor (SAS Institute 2003) 61 . Field trail data was transformed by ln (x + 1). Tukey's multiple-range test (α = 0.05) was used to identify the difference between cultivars in all the experiments. PCA was performed in Multibase 2013 in Microsoft Excel (www.numericaldynamics.com) by examining the correlation similarities between the observed measurements and four biological parameters and annual field population density of N. lugens were used as factors.  Table 3. Rice cultivars used for laboratory and field tests. BPH*, brown planthopper, N. lugens.