Introduction

Spiders are effective predators in terrestrial ecosystems, where they act in the ecology of herbivore and detritivore food webs (Jung et al. 2008a) and influence herbivore populations (Marc and Ysnel 1999). Spider species richness also is used in ecological conservation and assessment of environmental changes (Duelli 1998; Sauberer et al. 2004). The common wolf spider Pardosa pseudoannulata (Bosenberg et Strand) (Araneae: Lycosidae) is a major predator in Chinese agroecosystems (Tian et al. 2012). It is a biological control agent of some insect pests, especially the brown planthopper (BPH), Nilaparvata lugens (Hemiptera: Delphacidae), a severe rice pest in China (Zhao and Chen 2004; Zhao and Zhang 2005).

Rice, Oryza sativa, is a major staple food for more than half of the world’s population (Ouyang et al. 2005). However, the rice industry has been increasingly damaged by rice pests in the past few years (Bottrell 2012). Various traditional and biotechnology approaches have been employed to develop durable resistance to diseases and serious pests (Khush 2005). Genetically modified (GM) rice expressing the Bacillus thuringiensis (Bt) Cry1Ab protein has been developed to alleviate the environmental contamination due to the use of broad-spectrum chemical insecticides and to address the broader complex of insect pests (Chen et al. 2011).

Because the biotechnology of creating and deploying Bt rice is available, studies were launched to thoroughly investigate the agroecological and environmental safety of Bt rice. Akhtar et al. (2010) evaluated the non-target impact of Bt rice on the thrips, Stenchaetothrips biformis, showing longer larval and pupal development and previposition durations, but shorter oviposition period and female adult longevity and fewer total laid eggs of the thrip when fed on Bt rice in comparison to non-Bt controls (Akhtar et al. 2010). Chen et al. (2009) indicated that Bt rice did not influence survivorship and fecundity of Pirata subpiraticus but contribute to a longer developmental time of spiderlings (Chen et al. 2009). Lee et al (2014) indicated that the transgenic Cry1Ac rice had no adverse effects on the spider community structure of the rice fields (Lee et al. 2014). Lu et al (2014) assess the potential adverse effects of Bt rice on the non-target herbivore, Nephotettix cincticeps, and found there were no significant differences in some ecological fitness parameters including egg development duration, fresh adult weight, adult longevity, and oviposition period (Lu et al. 2014). However, the influence of Bt rice on behavioral ecology and its underlying molecular mechanisms in non-target organisms remain unexplored. In the present study, we performed transcriptomic analysis of CNSs isolated from control spiders maintained on non-transgenic rice and from experimental spiders taken from Bt rice, and identified differentially expressed genes (DEGs). We also recorded the influence of Bt on foraging behavior of spider under laboratory conditions, and found Bt rice caused a more active foraging behavior of the spider.

Materials and methods

Rice sample preparation

The seeds of the transgenic Bt Shanyou 63 rice, expressing a Cry1Ab gene, and its non-transgenic parental commercial cultivar Shanyou 63 were provided by Life Science College, Hunan Normal University. Both rice lines were grown in the experimental farmland at Hunan Agricultural University (113.08 °E, 28.18 °N). They were individually potted and cultured under the nylon nets (175 × 185 × 225 cm) to preclude extrinsic interference, such as unrelated insects and predators from the neighboring fields. No insecticides were sprayed in or near the field during the experimental period.

Animal sample preparation

Laboratory brown planthopper (BPH) colonies were established from the experimental field of Hunan Agricultural University. Each colony was reared on either Cry1Ab expressing (experimental) or control rice. BPHs were separately taken from both rice lines after 15 days continuous feeding, then provided to spiders as prey (Tian et al. 2013).

Adult female spiders (P. pseudoannulata) were collected from the organic farmlands of Hunan Academy of Agriculture Sciences and each spider was kept in a separate glass tube (15 mm diameter, 100 mm height). Each tube was topped with cotton balls and the inner tube wall was covered with a piece of soaked cotton to maintain humidity. Collected spiders were assorted into two groups (experimental and control) and starved for 48 h before being supplied each with 20 BHPs from experimental or control rice per day. All the tubes were numbered and remained in a temperature-controlled chamber (25 ± 2 °C, 70 ± 10% RH and 16 L:8D photoperiod; the photophase was from 06:00 to 22:00). In addition, our previous study indicated that accumulated Cry1Ab level in spider peaked at 9th day after spider was fed with Bt rice-consuming BHPs (Wang et al. 2016), so spiders used in this research were all fed with Bt containing BPHs for 9 days. Spiders were then not given food for 48 h before testing.

Effect of Cry1Ab rice on the foraging behavior of P. pseudoannulata

Fifteen spider adults of similar size (five spiders per repetition) were selected from the test and control groups respectively, transferred to individual cylinders (36.8 cm height, 6.5 cm diameter). The bottom of each cylinder was covered with moist mud (3 cm height) at the farmland of Hunan Agriculture University and a rice stem (30 cm height) provided by Hunan Academy of Agriculture Sciences, was vertically placed in the mud to simulate the farmland environment. Each cylinder was provided with 200 non Bt-rice reared BPHs and sealed with triple layers of gauze. The experiment was conducted in the temperature-controlled chamber as described and the foraging behavior was recorded by video camera for 48 h. The following periods (3rd–4th hour; 28th–29th hour; 35th–36th hour) from the treated and control groups were selected for data collection.

The movement behavior of the spider was recorded and analyzed using the following method: we marked the position of the spider in the cylinder every 10 s through video replaying. And the linear distance between two continuous points was calculated by means of Pythagorean proposition, which represented the movement distance of spider per 10 s (Fig. 1). We divided the range of moving distance of the spider into five categories, including 0–3 cm, 3–6 cm, 6–9 cm, 9–12 cm, and >12 cm. Then the occurrence frequency of each moving range was calculated during 1 h period.

Fig. 1
figure 1

The foraging displacements of P. pseudoannulata formed in a three-dimensional coordinate system. A spider was captured at the position A (X1, Y1, Z1) and B (X2, Y2, Z2) every ten seconds. A unit of foraging displacement was recorded by calculating the two ordinates absolute distance

Full size image

Transcriptomic analysis

RNA isolation and sequencing

CNSs were isolated from experimental and control spiders at 4 °C, transferred to 1.5 ml Eppendorf tubes, frozen with liquid nitrogen and stored at −20 °C. RNA extraction and sequencing was performed by Oebiotech Enterprise, Shanghai. In brief, total RNA was extracted from central nervous system (CNS) of spider, using TRIzol (Invitrogen Corp, USA) according to the manufacturer’s instructions. The integrity and quality of total RNA were checked using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc, Rockland, DE, USA), limiting the standard to 1.8 ≤ OD260/OD280 ≤ 2.1 and further confirmed by agarose gels.

RNA sequencing libraries were constructed and sequenced on an Illumina Hiseq 2000. Clean reads were assembled using the de novo transcriptome assembler Trinity after removing adaptor sequences, and low-quantity reads (reads with ambiguous bases “N” and duplicate sequences). Trinity combines relatively short, clean reads into longer continuous sequences, contigs, and assembles the contigs into longer unigenes. The unigenes were processed to remove redundancy, leading to non-redundant unigenes, which were used in the bioinformatics analysis reported here.

Gene annotation and expression

All unigenes (>200 bp) were analyzed against public databases (Nr/NCBI, SwissProt/UniProtKB, KOG/NCBI, KEGG, GO) using BLAST algorithm to understand gene biological functions. Three arthropod species databases (Tetranychus urticae, Drosophila melanogaster, and Stegodyphus mimosarum) were interrogated to assess identity of the obtained unigenes (e-values < 1e-5).

Expression of all unigenes in the spiders fed with test and control BPHs were estimated by calculating read density as ‘reads per kilobase of exon model per million map reads’ (RPKM) (Mortazavi et al. 2008). Employing the method by Audic and Claverie (1997), a strict algorithm was developed to identify differentially expressed genes (DEGs) between two libraries. A threshold false discovery rate of <0.01 and absolute value of log2 Foldchange >2 were used to determine significant differences in gene expression.

qPCR

To verify annotations of a set of selected genes, total RNA was isolated from each sample using TRIzol (Invitrogen Corp, USA) and treated with DNase I (Fermentas, Lithuania) according to the manufacturer’s protocol. The cDNA was synthesized using a RevertAid™ H Minus First Strand cDNA Synthesis Kit (Fermentas, Lithuania) and qPCR was performed using ABI 7900HT (ABI, USA) with a reaction volume of 25 μl, containing 1 μl of 1:10 cDNA diluted with ddH2O, 12.5 μl 2x SYBR Green Master Mix (ABI, USA) and 200 nM of primers. Reactions were performed in triplicate to ensure consistent technical replication and run in 96-well plates under the following conditions: 94 °C for 3 min. and then followed by 40 cycles for 94 °C for 30 s, 59 °C for 30 s, 72 °C for 45 s. Cycling ended at 72 °C for 5 min. Relative gene expression was evaluated with a DataAssist Software version 3.0 (Applied Biosystems/Life Technologies), using 18 s rRNA as a reference gene for RNA load and gene expression in analyses. The relative quantitative method (△△CT) was used to calculate the fold change of target genes (Livak and Schmittgen 2001).

Statistical analysis

The data from foraging behavior and qPCR were analyzed using one-way analysis of variance (ANOVA) with SPSS 17.0 software. Significant differences at p < 0.05 or p < 0.01 are designated with * or **.

Results

Effects of Cry1Ab protein on the foraging behavior of P. pseudoannulata

Control and experimental spiders were provided abundant control BPHs (n = 200). During an hour time-period, all recorded displacements were categorized into five distance zones (0–3, 3–6, 6–9, 9–12, >12 cm). The foraging range of P. pseudoannulata on a rice stem, as a main parameter, reflected their foraging capability. Spiders foraged most frequently at the 3–6 cm range, both in test and control groups, followed by 0–3 cm. Significant differences in the foraging range greater than 9 cm (p < 0.01) were observed between the test and control spiders (Fig. 2).

Fig. 2
figure 2

The total foraging displacements of P. pseudoannulata on the different movement range of rice within 1 h. The horizontal axis represents range of moving distance of spider per 10 s. The vertical line represents the occurrence frequency of each moving range during an hour. Mean (±SE) followed by asterisk (*) (**) in the figure were significantly different (One-way ANOVA, p < 0.05, p < 0.01)

Full size image

De novo assembly of all unigenes

The spider CNS cDNA sample was prepared and sequenced using the Illumina Hiseq 2000 platform. We obtained nearly 95 million total raw reads, which was submitted to the SRA database (accession number: SRR 2010603, SRR 2010603). Finally, a total of 120, 985 unigenes were obtained, using Trinity software (Table 1). The mean unigene size is 529.73 bp, with lengths ranging from 201 to 34055 bp. The size distribution of these unigenes is shown in Figure S1.

Table 1 Summary for the P. pseudoannulata CNS transcriptome
Full size table

Functional annotation

To date there remains no published complete reference genome for P. pseudoannulata. We put our sequences of CNS into five public databases and three arthropod databases as described in material and method section. Among all 120, 985 unigenes, 34,081 unigenes (28.17% of the total) were matched into Nr database, 29,472 unigenes (24.36% of the total) to Swissprot database, 27,234 unigenes (22.51% of the total) to COG, 29,194 unigenes (24.13% of the total) to GO database and 5,493 unigenes (4.54% of the total) to KEGG database. And 7,042 unigenes (5.82% of the total) showed similar sequences with the genome of Tetranychus urticae, 9,829 unigenes (8.12% of the total) with the genome of Drosophila melanogaster and 9,323 unigenes (7.71% of the total) with the genome of Stegodyphus mimosarum.

Overall, 27,234 unigenes were functionally classified into 25 COG categories (Fig. 3). Among them, the largest category was “Signal transduction mechanisms” (9864; 36.22%) followed “General function prediction only” (9399; 34.51%) and “Posttranslational modification, protein turnover, chaperones” (4581; 16.82%), while the categories with the least representation were “Cell motility” (125; 0.46%) and “Nuclear structure” (249; 0.91%).

Fig. 3
figure 3

COG classification for all unigenes

Full size image

The GO project provides an international standardized gene functional classification system for genes and gene products across all species. The ontology covers three domains: cellular component, molecular function and biological process (Ashburner et al. 2000). The unigenes assigned to the GO database were divided into three gene ontology classes: molecular function, cellular component and biological process. Among them, biological process made up the majority (64.49%), followed by molecular function (25.5%) and cell component (10.01%). The largest categories in gene ontology classes, “biological process”, “molecular function”, and “cell component” were “integral component of membrane”, “zinc ion binding”, and “oxidation-reduction process” respectively (Fig. 4).

Fig. 4
figure 4

GO classification for all unigenes. The Blast2GO program was used to obtain GO annotation of all unigenes (level 2). The top 10 GO terms which annotated most from three gene ontology classes, “biological process”, “molecular function”, and “cell component” were listed respectively

Full size image

We conducted a blast search against the KEGG database for the assembled unigenes. Results showed that 5,493 unigenes had significant matches in the database and were assigned into 327 KEGG pathways (Table S1). The pathway most represented by the unique sequences was “Ribosome”, “RNA transport” and “protein processing in endoplasmic reticulum”.

Functional analysis of DEGs

In all, the expression of 42 unigenes was significantly changed in Cry1Ab-treated spiders, compared to controls (FDR < 0.01, Log ratio >2, Table S2). Forty-two GO terms were annotated by DEGs, which could be divided into three parts including energy, binding and others (Fig. 5a). In specific, 14 GO categories were acting in regulating energy metabolism in intracellular, such as “oxidative phosphorylation” (GO:0006119), “cytochrome-c oxidase activity” (GO:0004129), “mitochondrial electron transport” (GO:0006120), mitochondrial membrane”(GO:0031966) (Fig. 5b); 8 GO categories were related to binding, such as “calcium ion binding” (GO:0005509), “iron ion binding” (GO:0005506), “heme binding” (GO:0020037), “copper ion binding” (GO:0005507), and “zinc ion binding” (GO:0008270) (Fig. 5c); 21 GO categories were involved in the regulation of diverse cellular processes, such as “proteolysis” (GO:0006508), “transport” (GO:0006810), “cysteine-type endopeptidase activity” (GO:0004197) and “peptidase activity” (GO:0008233) (Fig. 5d).

Fig. 5
figure 5

GO classification for DEGs. A  indicates the distribution of total GO terms, B indicates the distribution of terms of regulating energy metabolism, C indicates the distribution of terms of binding, D indicates the distribution of terms of regulation of diverse cellular processes

Full size image

COG classification also presented some categories related to energy metabolism, including “NADH dehydrogenase, subunit 4” (KOG4845),“NADH dehydrogenase subunits 2, 5, and related proteins (KOG4668)”, “Cytochrome c oxidase, subunit I” (KOG4769), “Cytochrome c oxidase, subunit II, and related proteins” (KOG4767) and “Cytochrome oxidase subunit III and related proteins” (KOG4664) (Table 2).

Table 2 COG classification for DEGs
Full size table

qPCR analysis

Four annotated unigenes (FDR < 0.05) that may have connection with the differences in the wolf spider foraging behavior were selected for qPCR analysis using primers listed in Table S3.

Uunigene comp16097_c0_seq1, similar to NADH dehydrogenase subunit IV, a key enzyme involved with the respiratory chain, was up-regulated by 1.5-fold. Unigenes comp76553_c0_seq1, comp50268_c0_seq1 and comp50316_c0_seq1, all akin to cytochrome oxidase subunit I, an enzyme related with electron transport chain, were up-regulated by 5.1-fold, 1.6-fold, 1.4-fold, respectively (Fig. 6).

Fig. 6
figure 6

qPCR analyses of the expression profiles of four annotated DEGs. Three technical replicates were performed for each of three biological replicates. The height of each box represents the mean average of sample-specific 2-DCt values. Bars displayed mean + SD and asterisks (*) (**) in the table are significantly different (One-way ANOVA, p < 0.05, p < 0.01)

Full size image

Discussion

Our previous study showed the multi-trophic movement of Cry proteins can be traced from Bt rice to BPHs and to a predatory spider (Wang et al. 2016), implying a potential negative influence on safety of predators in the field. In this research, we probed the impacts of Bt protein on behavioral ecology like foraging behaviors. Our foraging trial data show that Cry intoxication influenced spider foraging displacement where the movement range was greater than 9 cm. This implied that Cry1Ab intoxication has the potential to influence the foraging capability of P. pseudoannulata. In Apis mellifera, Cry1Ab affected food consumption and learning performance, hereby influencing honey bee foraging efficiency (Ramirez-Romero et al. 2008). Based on laboratory experiments with beet armyworms, Spodoptera exugua, Jiang et al. reported that low doses of Cry1Ac protein in an artificial larval diet led to increased long-flight behavior of beet armyworms (Jiang et al. 2013).

Usually, the CNS was a significant regulation core in organisms. In this study, we first conducted CNS transcriptome sequencing of P. pseudoannulata, investigating the molecular response of spider to Cry1Ab. Functional analysis of DEGs indicated that some genes were predicted regulating energy production in mitochondria, such as “oxidative phosphorylation”, “cytochrome-c oxidase activity”, “mitochondrial electron transport”, and “mitochondrial membrane”. These results are similar to that in Tenebrio molitor. The genes expression level involved in mitochondrial respiration was significantly increased in Tenebrio molitor larvae when they intake Cry3Aa protein (Oppert et al. 2012). The mitochondria is a powerhouse in cell, providing energy for various intracellular metabolic reaction (Conley 2016). Obviously, Bt proteins presented a negative impact on energy metabolism in cells. We speculated the changed energy metabolism in the spider may contribute to a more active foraging behavior than controls. On the other hand, some DEGs also have functions related to metal ion binding. Studies suggested that proteins (metalloproteins) contain a metal ion while metals are als required for proteins to carry out their functions (Waldron et al. 2009). These proteins are always acting in protein transportation, enzymes activity and signal transduction (Yi et al. 2009). For example, the iron is a common electron-transfer vector involved in the mitochondrial electron transport chain (Schuler et al. 1999). Hence, some metal ion-dependent reactions may be affected in response to Cry1Ab in CNS of P. pseudoannulata. This may be associated with a Cry1Ab tolerance mechanism in the spider.