Introduction Aggregation of animals is a common phenomenon that can be found in herds, shoals, flocks, schools and swarms of various species [1]. During aggregating, many species exhibit behavioral plasticity as an adaptive response to environmental changes [2]. Recently, several genes associated with behavioral plasticity have been identified, such as Cyp6a20 gene in fruit flies [3], foraging gene in honey bees [4], and Zenk gene in zebra finches [5], but the temporal dynamics of the interactions between genetic and environmental factors on behavioral plasticity are still largely unknown. The migratory locust, Locusta migratoria, is a worldwide agricultural pest, which exhibits severe plagues during which they migrate in marching juvenile bands and adult migratory swarms [6]. In response to the population density changes, they exhibit extreme phenotypic plasticity between the quiet, cryptically colored “solitarious” phase and the swarm-forming, conspicuously colored “gregarious” phase. Of all phase-related characteristics, behavioral changes are the most significant traits involved in the swarm formation in all locust species, particularly the changes in activity levels and in responses to other individuals [7]. Although the phase transition theory was originally established in the migratory locust by Uvarov [8], the best characterized model is the desert locust, Schistocerca gregaria. It is revealed that the behavioral gregarization in the desert locust can be triggered by a combination of sensory inputs from head (visual and olfactory) or by the mechanosensory stimulation of hind legs [9], [10]. These cues, which are stimulated by social environmental factors, are integrated by the nervous system to initiate behavioral phase change [11]. Recent genomic approaches have significantly improved our understanding in the phenotypic plasticity of the migratory locust. The application of large-scale expressed sequence tags (ESTs) revealed that 532 genes show different expression levels between the two phases of the migratory locust [12], [13]. Several hsp genes display higher expression levels in gregarious locusts and are regulated during phase change [14]. The two extreme phases exhibit different transcriptomic profiles of small RNAs, among which longer small RNAs, possibly piRNAs, are more abundant in the solitarious phase of the migratory locust [15]. As important sources of small RNA, transposable elements also show different expression patterns between the two phases, especially in central nervous and peripheral olfactory tissues [16]. However, to further understand the molecular mechanism of phase change, it is required to address the relationship between gene expression and behavioral changes, especially the time window when sensory cues and crowding trigger the behavioral transition. In this study, we developed a large-scale oligonucleotide microarray in the migratory locust and investigated the genome-wide expression profiles in heads of fourth-instar nymphs during the time courses of behavioral solitarization and gregarization. A binary logistic regression model was used to quantify the state of behavioral phase and the time course of behavioral phase change [17]. Furthermore, we used RNAi approach to explore the roles of candidate genes in phase-related behavior. The data from genome-wide expression profiling, the knockdown of olfactory-related genes and behavioral assays revealed the importance of olfactory-related genes, several CSPs and one takeout, in the fundamental phase-related behaviors of the migratory locust. Results The time course of behavioral phase change Based on the data collected from 86 gregarious and 69 solitarious fourth-instar migratory locust nymphs, we built a binary logistic model that can accurately categorize 87.2% of gregarious nymphs and 87.0% of solitarious nymphs, to summarize the behavioral phase state as a single probabilistic metric of gregariousness—Pgreg. The model retained three variables: attraction index, total distance moved and total duration of movement (Table 1). Using this model, we investigated the behavioral phase changes of fourth-instar nymphs during the time courses of solitarization (Isolation of gregarious locusts, IG) and gregarization (Crowding of solitarious locusts, CS). Gregarious nymphs showed significant behavioral solitarization after IG for only 1 h, which was in strong contrast to gregarious controls (Mann-Whitney U = 240, p = 0.001. Figure 1A). Full behavioral solitarization was acquired after IG for 16 h, exhibiting no significant difference with solitarious controls (Mann-Whitney U = 366.5, p = 0.110). However, solitarious nymphs did not display significant behavioral changes in gregarization until CS for 32 h as compared with solitarious controls (Mann-Whitney U = 261, p = 0.002. Figure 1B). Full behavioral gregarization was not acquired even after CS for 64 h as compared with gregarious controls (Mann-Whitney U = 115, p = 0.000). 10.1371/journal.pgen.1001291.g001 Figure 1 The time course of behavioral phase change in fourth-instar nymphs of Locusta migratoria. (A) Isolation of gregarious locusts (IG). (B) Crowding of solitarious locusts (CS). Pgreg, probabilistic metric of gregariousness. Arrows indicate median Pgreg values. 10.1371/journal.pgen.1001291.t001 Table 1 Variables retained in binary logistic model. Model Model if term removed Variables in equation β Standard error Wald df Sig. Exp(β) Model Log Likelihood Change in −2 Log Likelihood Sig. of the Change X1 attraction index 0.005 0.001 12.755 1 0.000 1.005 −56.557 15.875 0.000 X2 total distance moved 0.012 0.005 5.972 1 0.015 1.012 −52.005 6.770 0.009 X3 total duration of movement 0.015 0.007 4.694 1 0.030 1.015 −51.362 5.484 0.019 Constant −2.110 0.425 24.615 1 0.000 0.121 Pgreg = eη/(1+eη), η = −2.110+0.005×attraction index+0.012×total distance moved+0.015×total duration of movement. In addition, we also examined the time course of three behavioral variables that were significant predictors of phase-related behavior as revealed by logistic regression analysis. During the IG process, attraction index declined rapidly in 1 h to a stable level, but not to the level of solitarious controls even after IG for 64 h (Mann-Whitney U = 336, p = 0.037. Figure S1A). However, total distance moved and total duration of movement decreased to the level of solitarious controls in 4 h (Mann-Whitney U = 446, 375, p = 0.632, 0.092. Figure S1B and S1C). During the CS process, these three behavioral variables did not increase until 32 h (Mann-Whitney U = 318, 249, 292, p = 0.021, 0.001, 0.004. Figure S1D, S1E, S1F), but were far below the level of gregarious controls (Mann-Whitney U = 319, 220, 239, p = 0.035, 0.0004, 0.001, respectively). After CS for 64 h, only attraction index increased to the level of gregarious controls (Mann-Whitney U = 387.5, p = 0.359. Figure S1D). Gene expression profiling during the time course of phase change To investigate the global gene expression patterns related to behavioral phase change, we performed comparative gene expression profiling in heads of gregarious and solitarious nymphs that underwent population density changes. During the time courses of IG and CS, 794 and 1103 genes, respectively, were differentially regulated (Table S1, Dataset S1), and they fell into a diversity of ontological categories (Figure S2) and pathways (Table S2). The number of differentially expressed genes exhibited a wave-like pattern during the IG process, with the highest level at 16 h. The number of differentially expressed genes increased along with the time course of CS, with the highest level at 32 h (Figure 2A). Further hierarchical cluster analysis for all differentially expressed genes revealed fluctuating patterns for most genes during the IG process, but a relatively stable trend during the CS process (Figure S3). 10.1371/journal.pgen.1001291.g002 Figure 2 Candidate genes involved in locust phase change. (A) The numbers of up- and down-regulated genes during the time course of IG or CS. Numbers following IG or CS indicate hours of treatment. (B) Venn diagram indicating 453 mutually differentially expressed genes during the time courses of IG and CS. These 453 genes were categorized into 8 clusters by hierarchical clustering using a complete linkage algorithm. Cluster A and cluster H exhibit consistent changes over the time course of treatments and display a negative correlation of expression trends (Pearson correlation, r = −0.896, p 1.5 were considered as differentially expressed. Unsupervised hierarchical clustering [48] was performed with cluster 3.0 software using uncentered Pearson correlation and complete linkage, and presented by Java Treeview software [49]. Significant pathways were analyzed by KEGG Orthology-Based Annotation System (KOBAS, http://kobas.cbi.pku.edu.cn). InterPro categories were enriched for supplied gene list based on the algorithm presented by GOstat [50]. For each InterPro term, the difference between tested and reference gene groups was represented by a p value, which is estimated by chi square test. Fisher exact test was used when any expected value of count is <5. The p values of multitest were corrected by FDR [51]. Categories were considered as significantly enriched at p<0.05. Quantitative RT-PCR (qRT-PCR) cDNA was reverse-transcribed from 2 µg total RNA using MMLV reverse transcriptase (Promega). The standard curve method [52] was used to measure the mRNA relative expression levels, which were normalized by beta-actin [16], [42]. We adjusted all values of relative expression levels of these genes for all qRT-PCR experiments based on the starting concentration of plasmids initially used to dilute a series of standard curves in tissue-specific experiments. PCR amplification was conducted using Mx3000P spectrofluorometric thermal cycler (Stratagene) and RealMasterMix (SYBR Green) kit (Tiangen), initiated with a 2-min incubation at 95°C, followed by 40 cycles of 95°C, 20 s; 58°C, 20 s; 68°C, 20 s. Melting curve analysis was performed to confirm the specificity of amplification. RNA interference (RNAi) Double-stranded RNA (dsRNA) of green fluorescent protein (GFP), LmigCSP3 or LmigTO1 were prepared using T7 RiboMAX Express RNAi system (Promega) following the manufacturer's instructions. For both gregarious and solitarious locusts, fourth-instar locusts remained uninjected or were injected with 18 µg (6 µg/µl) of dsGFP, dsLmigCSP3, dsLmigTO1 or mixed dsLmigCSP3 and dsLmigTO1 (dsLmigCSP3& TO1) in the second ventral segment of the abdomen. Then, the injected gregarious locusts were marked and put back into gregarious-rearing cages, and the injected solitarious locusts were put back into solitarious-rearing cages. Three days later, the effects of RNAi on mRNA relative expression levels were investigated by qRT-PCR and the behavior was examined as described above. Statistical analysis Differences between treatments were compared either by Student's t-test or by one-way analysis of variance (ANOVA) followed by a Tukey's test for multiple comparisons. Behavior-related data were analyzed by Mann-Whitney U test because of its non-normal distribution characteristics. Chi square test was used for significance test between numbers of locusts in each arm in Y-tube assay. Differences were considered significant at p<0.05. Values were reported as mean ± SEM. Data were analyzed using SPSS 15.0 software (SPSS Inc., Chicago, IL). Supporting Information Dataset S1 Differentially expressed genes during time course of IG and CS. (0.47 MB XLS) Click here for additional data file. Figure S1 Change trends of three retained variables during the time course of IG or CS. Changes of attraction index (A), total distance moved (B) and total duration of movement (C) during the IG process. Changes of attraction index (D), total distance moved (E) and total duration of movement (F) during the CS process. Treatments are compared with gregarious control during the IG process or solitarious control during the CS process. *, p<0.05; **, p<0.01. (1.00 MB EPS) Click here for additional data file. Figure S2 Gene ontological categories of differentially expressed genes during the time course of IG or CS. Only differentially regulated gene categories (p<0.05) are shown here. (2.97 MB EPS) Click here for additional data file. Figure S3 Hierarchical clustering of differentially expressed genes during the time course of IG or CS. A complete linkage algorithm was used to cluster the 794 (IG) or 1103 (CS) differentially expressed genes. Horizontal stripes represent genes and columns show experimental treatments. Logarithmic fold change of treatment vs. reference (gregarious during IG or solitarious during CS) are shown in the heat map using red and green color codes for up- and down-regulation, respectively. (3.11 MB EPS) Click here for additional data file. Figure S4 Phylogenetic tree of representative CSP (A) and takeout (B) protein sequences. Consensus unrooted trees were generated with 1000 bootstrap trials using the neighbor-joining method and presented with a cutoff value of 50. (0.98 MB EPS) Click here for additional data file. Figure S5 Alignment of deduced protein sequences of CSPs from locust EST database. Putative signal peptides (shaded in yellow) were predicted using SignalP V3.0 web server, and conserved cysteine residues and sequence regions are indicated by red letters. (6.31 MB EPS) Click here for additional data file. Figure S6 Alignment of deduced protein sequences of takeout from locust EST database. Putative signal peptides (shaded in yellow) were predicted using SignalP V3.0 web server, and conserved cysteine residues are indicated by red letters. (7.91 MB EPS) Click here for additional data file. Figure S7 Tissue-specific expression of LmigCSP1, LmigCSP2, and LmigCSP4 in gregarious and solitarious fourth-instar nymphs. A, antenna; L, labial palp; B, brain; W, wing; HL, hind leg; FB, fat body. Means labeled with the same letter within each treatment are not significantly different and error bars represent SEM. *, p<0.05; **, p<0.01. (0.87 MB EPS) Click here for additional data file. Figure S8 Gene expression of three other CSP genes in antennal tissue during the time course of IG or CS. (A) RNA relative expression levels of LmigCSP1, LmigCSP2 and LmigCSP4 during the time course of IG (ANOVA, F5,18 = 2.810, 18.346, 9.439, p = 0.066, 0.000, 0.001, respectively). (B) RNA relative expression levels of LmigCSP1, LmigCSP2 and LmigCSP4 during the time course of CS (ANOVA, F5,18 = 1.474, 4.827, 8.020, p = 0.269, 0.012, 0.002, respectively). (0.82 MB EPS) Click here for additional data file. Figure S9 Systemic effects of RNAi. (A) RNA relative expression levels of LmigCSP3 (Upper) and LmigTO1 (Lower) genes after injection of dsLmigCSP3 and dsLmigTO1, respectively, in the four selected tissues in gregarious nymphs. (B) RNA relative expression levels of LmigCSP3 (Upper) and LmigTO1 (Lower) genes after injection of dsLmigCSP3 and dsLmigTO1, respectively, in the four selected tissues in solitarious nymphs. A, antenna; L, labial palp; B, brain; HL, hind leg. In all groups, treatment (inj-dsGFP, inj-dsLmigCSP3 or inj-dsLmigTO1) compared with respective non-injected controls. *, p<0.05; **, p<0.01; n.s., not significant. (0.99 MB EPS) Click here for additional data file. Figure S10 Six CSP genes have a high sequence identity of 78.18%. GenBank accession number: LmigCSP5, CO835786; LmigCSP6, CO852124. (2.12 MB EPS) Click here for additional data file. Figure S11 Effect of dsLmigTO1 injection on RNA relative expression levels of LmigTO2 and LmigTO3. Student's t-test, t = 0.372, 0.770, p = 0.723, 0.471, respectively. In all groups, treatment (inj-dsGFP or inj-dsLmigTO1) compared with respective non-injected controls. *, p<0.05; **, p<0.01; n.s., not significant. (0.66 MB EPS) Click here for additional data file. Figure S12 Interaction of LmigCSP3 and LmigTO1. (A) RNA relative expression level of LmigCSP3 after dsLmigTO1 injection in gregarious or solitarious nymphs. (B) RNA relative expression level of LmigTO1 after dsLmigCSP3 injection in gregarious or solitarious nymphs. In all groups, treatments (inj-dsGFP, inj-dsLmigCSP3 or inj-dsLmigTO1) are compared with respective non-injected controls. n.s., not significant. (0.81 MB EPS) Click here for additional data file. Figure S13 The behavioral phase state of gregarious nymphs (A) and solitarious nymphs crowding for 4 instars (B). Pgreg, probabilistic metric of gregariousness. Arrows indicate median Pgreg values. (0.63 MB EPS) Click here for additional data file. Figure S14 Schematic diagram of behavioral assay system. (A) Behavioral assay arena. (B) Y-tube olfactometer assay set. Individual nymph was released from Y-tube adapter to observe their choice. (8.20 MB EPS) Click here for additional data file. Figure S15 Schematic diagram of experimental design of behavioral assay and microarray hybridization. IG (left of the dashed), isolation of gregarious locusts; CS (right of the dashed), crowding of solitarious locusts. Numbers following IG or CS represent hours of treatment. Arrows indicate microarray hybridization of neighboring treatments. (1.01 MB EPS) Click here for additional data file. Table S1 Number of differentially expressed genes during the time course of IG or CS. (0.05 MB DOC) Click here for additional data file. Table S2 Top ten regulated pathways during the time course of IG or CS. (0.07 MB DOC) Click here for additional data file. Table S3 Genes used in phylogenetic tree reconstruction. (0.07 MB DOC) Click here for additional data file. Table S4 qRT-PCR validation of microarray data and primer sequences. (0.17 MB DOC) Click here for additional data file.
