A multi-generational risk assessment of Cry1F on the non-target soil organism Folsomia candida (Collembola) based on whole transcriptome profiling

Test organism

The Danish strain of the parthenogenetic F. candida Willem, 1902, was originally obtained from Aarhus University, Denmark, and has been cultured in our laboratory for over 10 years. F. candida is usually fed granulated dried baker’s yeast (Hebei Mauri Foods Co., Ltd., Zhangjiakou, China) and was reared in the Petri dishes (90 × 18 mm) containing a five to seven mm layer of a solidified mixture of plaster of Paris and activated charcoal (9:1 wt/wt, dissolved in distilled water, with about 270 ml of water/500 g of mixed powder) (Figs. 1A–1C). The baker’s yeast was placed on the surface of the plaster and was renewed weekly to reduce the growth of other fungi. Distilled water was added to the base as needed such that free water was present in the plaster pores but did not form a film on the plaster surface. These dishes with F. candida were kept in an artificial climate chamber in total darkness (20 ± 1 °C, 80% relative humidity).

Dietary exposure experiment

The artificial diet containing Cry1F (BT diet) was prepared as follows: two mg of purified Cry1F protein (Envirologix Inc., Portland, ME, USA) and four g of baker’s yeast granules were dissolved in 10 ml of distilled water; the preparation was fully mixed and then transferred into a plastic container. The final concentration of Cry1F protein in the BT diet was 500 μg/g, which is much higher than the Cry1F protein concentrations estimated in IRGE crops (less than 50 μg/g, ILSI Research Foundation, 2013), as well as than the EC₅₀ values of Cry1F against lepidopteran pests Chilo suppressalis (6 μg/g) and larvae of Bombyx mori (136 ng/g) (Jiao et al., 2016). After 27 h of lyophilization, the BT diet was ground into powder and stored at −70 °C. As the control, baker’s yeast granules without Cry1F (CT diet) were prepared in the same way. We also tried a positive control by adding the insecticidal compounds (potassium arsenate) into the baker’s yeast granules, but it significantly affected the survival and the development of F. candida, which hampered the performance of the consecutive multi-generational exposure experiment. Therefore, we didn’t perform the positive control in our multi-generational DEE. However, the concentration of Cry1F in the BT diet was high enough to validate the test of Cry1F to F. candida.

Based on the Organization for Economic Cooperation and Development test protocol (Organisation for Economic Cooperation and Development (OECD), 2016), we designed a multi-generational DEE on F. candida (Fig. 1D). During the experiment, all F. candida were cultured in the same artificial climate chamber with the culture conditions mentioned above and were synchronized to the same stage (neonates) at the start of each generation in the experiment.

First, we transferred 80–100 randomly selected F. candida adults (about 4–6 weeks old) to a new Petri dish. After 72 h of oviposition, all adults were removed. The eggs began to hatch 10 days later. After 24 h of hatching, we selected synchronized neonates. They were transferred to a new Petri dish and were fed with pure baker’s yeast for 10 days. Next, all juveniles were separated into two sets: the BT set were fed BT diet, and the CT set were fed CT diet. Each set was done with individually- and group-reared collembolans.

In the first-generation DEE, 60 juveniles were transferred to 60 small plastic Petri dishes (55 × 14 mm, plaster height: 5 mm), respectively, as individually-reared collembolans. A total of 30 were fed the BT diet, and 30 were fed the CT diet. The remaining juveniles were evenly distributed into six plastic Petri dishes (85 × 16 mm, plaster height: 5 mm), with 60–80 collembolans per dish. Three of them were fed the BT diet and the other three were fed the CT diet. All diets were renewed every 2 days to avoid Cry1F protein degradation and the growth of other fungi.

The first generation was assessed after 25 days of DEE, which is 35 days from hatching. We choose this time point to ensure that F. candida had sufficient exposure to Cry1F, and had laid one to two batches of eggs, but with only one batch of eggs hatching. For the individually-reared collembolans, the numbers of adults and juveniles in each treatment were counted to calculate the survival and reproduction rates. For the collembolans fed in groups, about 20 individuals per replicate were randomly collected for enzyme-linked immunosorbent assay (ELISA, see section Enzyme-linked immunosorbent assay below), and 25–35 individuals per replicate were collected for RNA-Seq (see section RNA extraction, sequencing, and unigene annotation below).

The remaining collembolans fed in groups of CT and BT sets were used for synchronization of the second generation. To maintain the same experimental conditions, the CT and BT sets were synchronized independently, and their neonates were fed the CT diet or BT diet directly after birth. After 10 days, the second-generation DEE began and was conducted in the same manner as the first-generation DEE except that there were 40 replicates per set for individually-reared collembolans and two replicates per set for group-reared collembolans. After 25 days of DEE, survival and reproduction were evaluated using the individually-reared collembolans. ELISA and RNA-Seq were performed using the group-reared collembolans as described for the first generation. The remaining collembolans that were fed in groups in the second generation were used to start the third generation. The procedure for the third-generation DEE was identical to that for the second except that there were three rather than two replicates per set for the group-reared collembolans.

Survival and reproduction

The individually reared collembolans were used to evaluate the survival and reproduction of F. candida at the end of each generation. Several replicates were excluded from the analysis: (1) missing adult collembolans (five in BT treatment and five in CT treatment over three generations), which are assumed to have escaped because there was no corpse left; (2) eggs didn’t hatch due to overgrown fungi (six in BT treatment and five in CT treatment over three generations).

With the aid of a stereoscopic microscope (SMZ-10; Nikon, Minato, Japan), the number of living adults was determined as a measure of survival, and the number of juveniles produced was determined as a measure of reproduction. We calculated the mean number of juveniles produced per individually fed collembolan. Due to the wide variation in the reproduction data (7–41 juveniles per individually-reared collembolan), the data were log-transformed prior to conducting statistical analysis. Student’s t-test was used to analyze the differences between the two treatments in each generation.

In addition, the data for the reproduction of all BT and CT sets over three generations were subjected to a two-way ANOVA, followed by the least significant difference test, with diet treatment, “generation,” and their interaction as fixed factors. Generation is not an independent variable, but we used “generation” as a label representing unknown variable factors in order to test whether the reproductions of F. candida over three generations in the same treatment were affected by unpredictable environment stress. All of these statistical analyses were performed with IBM SPSS Statistics 24 (version R24.0.0.0). Differences were considered significant at p < 0.05.

Enzyme-linked immunosorbent assay

The concentration of Cry1F in the diets and in the F. candida fed in groups was measured by ELISA for each replicate in each generation. The QuantiPlate Kit for Cry1F (Envirologix Inc., Portland, ME, USA) was used to detect Cry1F in a two- to four- mg sample of the fresh diet in each replicate, in a two- to four- mg sample of the diet in each replicate after 2 days of feeding, and in 20 individuals of F. candida (2–5 mg) per replicate. Before the test, all F. candida samples were washed in a PBST extraction buffer (phosphate-buffered saline with Tween-20, pH 7.4) to remove any Cry1F toxin from the outer surface of the collembolan cuticles. Then, all diet samples and F. candida samples were fully ground with an electric grinding rod and extracted with a PBST extraction buffer (phosphate-buffered saline with Tween-20, pH 7.4), respectively. ELISA was performed according to the manufacturer’s instructions. A two-way ANOVA was used to compare Cry1F concentrations of three BT sets over three generations by using the software IBM SPSS Statistics 24 (version R24.0.0.0). Differences were considered significant at p < 0.05.

RNA extraction, sequencing, and unigene annotation

The miRNeasy® Mini Kit (Qiagen, Hilden, Germany) was used to extract the total RNA from 25–35 individuals per replicate of the collembolans that were fed in groups for each generation. The RNA concentration and integrity was evaluated with Aglient 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).

The cDNA library construction, Illumina sequencing, and de novo assembly of RNA samples were carried out by Hangzhou 1gene Technology Co., Ltd (Hangzhou, China). All cDNA libraries were constructed using the NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) following the protocol described by the manufacturer. The libraries were sequenced with 150-bp paired-end reads on an Illumina Hiseq 4000 platform (San Diego, CA, USA). After sequencing, the raw reads (about 6 GB of data for each replicate) were filtered to remove adaptor sequences, duplication sequences, and low-quality sequences. The clean reads were de novo assembled into unigenes using Trinity and SOAPdenovo-Trans. Six public databases were used to annotate unigenes with BLASTx or BLASTn (E-value < 10⁻⁵), including NCBI Nr and Nt databases (http://www.ncbi.nlm.nih.gov/), SwissProt (http://www.expasy.ch/sprot/), KEGG (http://www.genome.jp/kegg/), COG (http://www.ncbi.nlm.nih.gov/COG/), and GO (http://www.geneontology.org/).

Differential gene expression and biological marker expression

Based on strict synchronization and multiple repeats of samples, we screened DEGs between groups of CT and BT sets in the same generation. The expression levels of all unigenes were calculated using the fragments per kb per million fragments (FPKM) method: FPKM = (10⁶C)/(NL/10³). The data of two or three replicates for each set in every generation were averaged, and the whole transcriptome gene expression of all CT and BT sets were compared. DEGs were determined by the condition of (fold change, |log2 ratio| of FPKM (BT vs. CT)) ≥ 1 and false discovery rate ≤ 0.05. Co-expressed DEGs over three generations were screened and displayed with Venn diagrams. All DEGs over three generations were analyzed by heatmaps and the hierarchical clustering.

Additionally, we carefully checked the expression of seven known biological markers, whose changes in expression level are often used to indicate that the test organism is experiencing stress. These genes encode the following proteins: the antioxidant-related enzymes catalase (CAT) (Niu et al., 2017; Yousef, Abdelfattah & Augustyniak, 2017) and superoxide dismutase (SOD) (Yang et al., 2015; Yousef, Abdelfattah & Augustyniak, 2017), the detoxification-related enzymes glutathione S-transferase (GST) (Niu et al., 2017; Oliveira et al., 2015), carboxylesterase (CES) (Niu et al., 2017; Yang et al., 2015), and glutathione reductase (GR) (Yang et al., 2015), a metallothionein-like motif containing protein (MTC) (Nakamori & Kaneko, 2013), which is unique in F. candida and which is very sensitive to different heavy metal ions, and the heat shock protein 70 (HSP70) (Liu et al., 2010), which is sensitive to insecticides, drought, and other environmental stresses. The expression profiles of unigenes annotated to these biological markers were examined and their expression levels were compared among different samples. All these genes of the seven biological markers, over three generations, were analyzed by heatmaps and the hierarchical clustering.

The software R (R Core Team, 2017) was used for the heatmaps (R package “pheatmap” (Kolde, 2019)) and the hierarchical clustering (“dist” and “hclust” orders, distance method as “euclidean,” and cluster method as “complete”) of all DEGs and biological markers over three generations (expression data of DEGs were transformed by using the common logarithm of the counts plus 1).

Cry1F concentrations in the diet and F. candida

Enzyme-linked immunosorbent assay measures showed that no Cry1F toxin was detected in the pure yeast diet or in the F. candida that were fed the CT diet. The average concentration (mean ± SE) of Cry1F was 415 ± 15 μg/g in the fresh BT diet and 358 ± 57 μg/g in the BT diet after 2 days of feeding, indicating that Cry1F was continually present in the experiment. The average concentrations of Cry1F in collembolans that were fed the BT diet were 3.5 ± 1.3, 2.7 ± 1.4, and 5.2 ± 1.0 μg/g at the end of the first, second, and third generation, respectively. These concentrations were not significantly different based on a two-way ANOVA (p = 0.429) and suggest that the Cry1F was continuously ingested by collembolans.

Survival and reproduction of F. candida

All collembolans in both the CT and BT sets survived and reached the adult stage for all three generations, suggesting that Cry1F has no obvious negative effect on collembolan survival. Over the three generations, the mean (± SE) number of juveniles produced per individually fed collembolan was 22.57 ± 0.635 for 100 samples in the CT set and 20.70 ± 0.590 for 99 samples in the BT set. According to Student’s t-test for each generation, the number of juveniles produced per collembolan did not significantly differ between the CT and BT sets for any of the three generations (p = 0.912, 0.114, 0.071 for generations 1, 2 and 3, respectively) (Fig. 2). Two-way ANOVA (with diet treatment, “generation,” and their interaction as fixed factors, Table 1) showed that the mean number of juveniles produced per collembolan was significantly affected by “generation” (p = 0.001) but not by diet treatment (p = 0.068) or the interaction between diet treatment and “generation” (p = 0.410). The results suggest that even in the same treatment culture condition, the reproductions of F. candida were still affected by some unpredictable environment stress instead of the BT diet.

Transcriptome sequencing and annotation

A total of 316,693,647 raw reads for 16 samples were generated, and deposited in the NCBI Sequence Read Archive (SRP132745). After data filtering and de novo assembly, 284,174,422 clean reads were assembled into 93,976 unigenes. The total length of these unigenes was 155,829,628 bp. The mean length of the unigenes was 1,658 bp, and N50 length was 3,292 bp. The sequence data of each type of sample were unbiased, and enough reads were obtained to perform gene expression analyses. The quantity and quality of the RNA-Seq data are shown in Table 2.

Of the total number of unigenes, 55,390 (58.94%), 17,230 (18.33%), 47,734 (50.79%), 43,284 (46.06%), 31,016 (33.00%), and 18,603 (19.80%) were annotated in NCBI Nr, NCBI Nt, SwissProt, KEGG, COG, and GO databases, respectively. Overall, 57,758 (61.46%) unigenes were annotated to known protein/nucleotide sequences.

Differentially expressed genes (DEGs) over three generations

A total of 463 DEGs (0.49% of all unigenes) between the CT and BT sets were identified for all three generations (Figs. 3A–3C), including 211, 19, and 244 DEGs for the first, second, and third generation, respectively (Fig. 3D). There was no consistent tendency for either the up-regulation or down-regulation of DEGs and there was no consistent co-expression of DEGs over the three generations (Fig. 3E) and only six DEGs (five up-regulated and one down-regulated) were identified in both the first and third generations. Five contra-regulated DEGs were detected in the second and third generations. The functions of these 11 DEGs are unknown.

The heatmap analysis of DEGs over three generations showed similar expression patterns of the CT and BT sets in each generation (Fig. 4A). The hierarchical clustering of all samples demonstrated that the CT and BT samples of the same generation were always clustered together, instead of the CT or the BT sets of different generations (Fig. 4B).

Biological marker expression profile

A total of 855 unigenes for the seven biomarkers (CAT, GST, SOD, GR, CES, MTC, and HSP70) were annotated, but only five unigenes were differentially expressed between the CT and BT sets, that is, were up- or down-regulated by more than two-fold in the BT set (Table 3): one unigene of CAT (CL2174.Contig5_All) was up-regulated in the third generation, three unigenes of CES (Unigene7933_All, CL5034.Contig2_All, and CL6365.Contig1_All) encoding CES type B were up- or down-regulated in the first or the third generations, one unigene of HSP70 (Unigene15285_All) was greatly down-regulated, that is, there was a five-fold decrease in expression, in the first generation. However, the expression of the five DEGs was significantly different in the BT vs. CT sets in only a specific generation, and none of those changes are consistent in all three generations (Table 3). Furthermore, the gene expression hierarchical clustering of the seven markers (Fig. 5) clearly showed that the CT and BT sets from the same generation always clustered together, which confirms that Cry1F treatment did not affect the functions related to the seven biological markers.