A Drosophila systems model of pentylenetetrazole induced locomotor plasticity responsive to antiepileptic drugs

Fly culture and harvesting

D. melanogaster wild type Oregon-R strain was used. Routine cultures were maintained at 24 + 1°C, 60% RH, and 12 hrs light (9 AM to 9 PM) and 12 hrs dark cycle. Standard fly medium consisting of agar-agar, maize powder, brown sugar, dried yeast and nipagin was used. Standard methods of fly handling and manipulation were followed. Stringency required in behavioral studies was strictly adhered to at various levels, conditions of housing, exposure to anesthetic agent, light intensity, for example. Three to four days old unmated male flies were used to begin all experiments except dsRNA microinjection. For microinjection, 5–6 day old flies were used.

Drug dosage and treatment

LC₅₀, defined as concentration in normal fly media causing 50% lethality in seven days, was determined for PTZ and the AEDs (results not presented). Unless mentioned otherwise, half of LC₅₀, 8, 3.48, 0.33 and 5 mg/ml for PTZ, ETH, NaVP (all from Sigma-Aldrich), and LEV (Levesam 500, Nicholas Piramal), respectively, was used. As GBP and VGB (both from Sigma-Aldrich) did not cause lethality up to 16 and 24 mg/ml respectively, the same concentrations were used in all the experiments. Thirty flies were housed in each treatment vial. Flies were maintained at 24 + 1°C, 60% RH, and 12 hrs light (9 AM to 9 PM) and 12 hrs dark cycle.

Locomotor assays

Climbing speed was measured using an indigenously developed semi-manual method that was validated by DIAS (v. 3.4.2, Soll Technologies). In semi-manual assay, a glass column of 2 cm internal diameter and 30 cm internal length, i.e., length between the two cotton plugs, was used in the assays. The column was marked with lines at every cm along the length. Each fly was first familiarized in the vertically placed column for 90 seconds. The fly was thereafter brought to the bottom of the tube, by tapping the tube on a piece of packing foam. As soon as the fly had fallen on the cotton plug at the bottom, the tube was as such placed vertically and, in semi-manual method, a "dot/comma" recording, explained below, was applied the moment fly crossed a height of 5 cm. In "dot/comma" recording, the locomotor activity of a single fly was monitored by keeping pressed the dot key or the comma key of a personal computer, to record a moving or resting fly respectively. Dots and commas were subsequently transformed in activity and rest period respectively, using the cursor speed. Climbing speed was calculated using the following formula, s = h/t, where s = climbing speed, h = height climbed in cm, and t = activity period in sec. One fly each from different treatment groups was first scored as described below, discarded, and the exercise repeated. Measurements were taken at room temperature, between 9 AM to 9 PM. Climbing was considered complete when the fly reached the cotton plug at the top, fell to the bottom after climbing a certain height beyond 5 cm mark, or stopped after climbing to a certain height for more than 10 sec. Comma was applied only when the fly took rest. The vertical locomotor assay was adapted to measure horizontal (spontaneous) locomotor activities. A single fly was first brought to the middle of the column by gentle shaking and then the fly movement was constantly monitored for 90 sec by keeping pressed the dot key or the comma key of a personal computer, to record a moving or resting fly respectively. The total number of lines that a fly crossed in the 90 sec assay time was counted and noted down at the end of the dot/comma recording. The data was normalized for 90 sec and the dots and commas were subsequently transformed in rest and activity period respectively, using the cursor speed. Walking speed was obtained by dividing the lines (cm) a fly crossed (distance walked, d) by time, in sec, it spent in activity (activity period, a), during the 90 sec assay period. For some flies, both "d" and "a" were found to be zero. Such individuals were excluded while calculating population mean of s. Finally, the horizontal assay was used to extract four locomotor parameters, namely, activity period, distance walked and walking speed. The semi-manual climbing assay described above was validated using DIAS. TD was measured using DIAS, by video recording flies in the treatment vials with media containing 16 mg/ml PTZ. TD is the net path length divided by the total path length. This gives 1.0 for a completely straight path and a smaller value for a meandering path.

Microarray expression profiling

Total cellular RNA was isolated from fly heads belonging to five biological replicates, four for microarray profiling, and one for real-time PCR. Each replicate represented four independent sets of control and treated flies. Each independent set consisted of eight vials, four each of NF and PTZ, treated in parallel. Similar design was followed for control NF-NF control microarray experiment. Thirty flies were treated in each vial. Heads of 120 flies were pooled for isolating each RNA sample. Flies were frozen in 50 ml falcon tubes in liquid nitrogen. Two cooled sieves were arranged such that the larger sieve (mesh size 850 mm) was placed on top of the smaller one (mesh size 355 mm). Frozen flies were shaken 4–5 times in the falcon and poured onto the top sieve. The flies were brushed gently with a paint brush till all heads were sieved out on the bottom sieve. Bodies that remained on the top sieve were discarded and heads were collected in cryovials and kept frozen at -80°C. Total RNA was isolated from frozen fly heads using TRI REAGENT (Sigma), according to the manufacturer's protocol. RNA was quantified using spectrophotometer. RNA quality was checked by running 1% agarose gel.

Microarray -cDNA Synthesis Kit, -Target Purification Kit, and -RNA Target Synthesis Kit (Roche) were used to generate labeled antisense RNA. Starting with 10 μg of total cellular RNA, Eberwine method (kits from Roche) was used to generate cDNA and thereafter Cy³ and Cy⁵ (Amersham) labeled antisense RNA. The Cy³ and Cy⁵ labeled aRNAs (control and treated) were pooled together and precipitated, washed, air-dried, and dissolved in 18 MΩ RNAase free water (Sigma). The recovery of labeled aRNAs was checked using spectrophotometer and agarose gel electrophoresis. A total of 16 microarrays (12Kv1, CDMC) were hybridized, four each for NF-NF control (0 hr); NF-PTZ, 12 hrs; NF-PTZ, 2^nd day; NF-PTZ, 7^th day. Out of four, two slides were dye-swaps. Hybridization solution was prepared by mixing hybridization buffer (DIG Easy Hyb, Roche), 10 mg/ml salmon testis DNA (0.05 mg/ml final concentration, Sigma) and 10 mg/ml yeast tRNA (0.05 mg/ml final concentration, Sigma) and added to the labeled product. This mixture was denatured at 65°C and applied onto cDNA microarray slides. The slides were covered by lowering down a 24 × 60 mm coverslip (ESCO, Portsmouth, USA). Hybridization was allowed to take place in hybridization chamber (Corning) at 37°C for 16 hrs. After hybridization, coverslips were removed by submerging the slides in a solution containing 1× SSC and 0.1% SDS at 50°C. Slides were washed (three times for 15 minutes each) in a coplin jar at 50°C with occasional swirling and then transferred to 1× SSC and washed with gentle swirling at room temperature (twice for 15 minutes each). Finally, slides were washed in 0.1× SSC for 15 minutes and then liquid was quickly removed from the slide surface by spinning at 600 rpm for 5 minutes. Slides were scanned at 10 μm resolution using GenePix 4000A Microarray Scanner (Molecular Devices), and the images preprocessed and quantified using Gene Pix Pro 6.0 (Molecular Devices).

Statistics

Unless mentioned otherwise, pair-wise Student's t-test, two-tailed, heteroscedastic, uncorrected, was performed for behavioral analysis. In microarray gene expression analysis, ratio based data normalization and selection of features were performed using Acuity 4.0 (Molecular Devices). All Spots with raw intensity less then 100 U and less then twice the average background was ignored during normalization. Normalized data was filtered for the selection of features before further analysis. Only those spot were selected which contained only a small percentage (< 3) of saturated pixels, were not flagged bad or found absent (flags > 0), had relatively uniform intensity and uniform background [Rgn R2 (635/532) > 0.6] and were detectable above background (SNR > 3). Analyzable spots in at least three of four biological replicates performed were retrieved for downstream analysis using SAM 3.0 (Excel Add-In) [38], under the conditions of one class response and 100 permutations. Wherever available, non-CG number genes were converted to CG numbers mainly using BDGP http://www.flybase.org. Gene IDs were converted to FBgn numbers using GeneMerge http://genemerge.bioteam.net/convertgenenames.html, before GOTool Box [39] was used to retrieve overrepresented biological processes in up- or down- regulated genes, under the settings, ontology, biological process; mode, all terms; reference, genome; evidence, all-all evidence; species, D. melanogaster; GO-stats; statistic test, hypergeometric http://burgundy.cmmt.ubc.ca/GOToolBox/. DAVID [40] was used for pathway enrichment analysis http://david.abcc.ncifcrf.gov/home.jsp. P values were obtained with or without correction for multiple hypotheses testing, as indicated in the results section. Hypergeometric distribution probabilities for genes and GO processes were calculated assuming population sizes of 10000, approximately the number of unique genes in the arrays, and 4041, a previously estimated number for Drosophila [41], respectively.

Chronic PTZ locomotor model and its validation by AEDs

We first explored, through empirical testing, if PTZ could produce proconvulsant-like effect in Drosophila adults. Though convulsions were not observed, visual examination indicated that flies treated with 16 mg/ml of PTZ for 12 hrs exhibit hyperkinetic behavior. Software-based analysis of video recordings showed a decreased 'total directionality' in PTZ treated flies compared to control, i.e., flies treated with normal food (NF). Mean + S.E of 'total directionality' was found to be 0.39 + 0.07 (n = 11) and 0.19 + 0.05 (n = 11) for flies fed with NF and with food containing PTZ, respectively. A significant (p = 0.04) difference in 'total directionality' demonstrated that 16 mg/ml PTZ causes flies to take a more circuitous path. As 16 mg/ml was the LC₅₀ dose (lethality data not presented), we used half of this dose, i.e., 8 mg/ml, as a sub-hyperkinetic dose for further developing and evaluating a chronic model. Though seizure-like behavior was never observed in flies treated with 8 mg/ml in a chronic manner, we empirically determined a regime in which seven days of chronic PTZ treatment and seven following days of PTZ discontinuation respectively caused a decrease and an increase in startle-induced climbing speed in flies (Figure 1; for DIAS validation of climbing assay, see additional file 1). Examining the effect of PTZ on various horizontal locomotor activities and non-locomotor parameters, namely, courtship duration, fertility, and body weight showed that chronic treatment induced decrease and subsequent withdrawal induced increase in climbing speed are specific (for detailed results, see additional file 2).

We next examined if AEDs are effective in ameliorating climbing alterations induced by chronic PTZ and withdrawal from chronic PTZ. Concomitant treatment with PTZ and ETH, GBP, VGB, NaVP or LEV for seven days showed the last two AEDs as effective in alleviating PTZ induced climbing speed deficit in flies on 7^th day (Figure 2A). NaVP and LEV were thus found to be effective in the fly model. To examine long term effectiveness of AEDs, we concomitantly treated flies with PTZ and an AED for seven days and measured climbing speed seven days after withdrawal of the two drugs. Only LEV, not other AEDs, was found to ameliorate climbing speed increase caused by withdrawal from chronic PTZ (Figure 2B). LEV thus exerted a long term effect. To examine the possibility that the above validation results might have been confounded by AEDs own locomotor effects, we treated flies with AEDs for seven days and measured climbing speed at the end of this chronic treatment as well as seven days after AEDs were withdrawn. ETH, VGB and LEV, not NaVP and GBP, were found to decrease climbing speed at the end of chronic treatment (Figure 3A). Independent speed lowering effect of NaVP and LEV (Figure 3A) could not be considered to confound the observation that these AEDs ameliorate decrease in climbing speed caused by PTZ (Figure 2A). After withdrawal, none of the AEDs were found to cause climbing abnormality (Figure 3B). Cumulatively, the above results validated the PTZ model as predictive of AEDs.

NaVP and LEV, not ETH, GBP and VGB, were found to be effective in our chronic PTZ model of locomotor plasticity. NaVP and LEV are known to suppress both kindling development as well as kindled seizure in rodent models [42, 43]. Clinical trials however suggest that LEV, not NaVP, possesses antiepileptogenic activity [9, 44, 45]. Overall, the profiles of NaVP and LEV in the fly model were consistent with their therapeutic potential. Ineffectiveness of other AEDs in the fly model however suggested that our model is only partly predictive of therapeutic potential. Clinical relevance of the model is also limited by the fact that the drug doses used in developing the model were selected either empirically or based on toxicity.

Time series of transcriptomic changes underlying chronic PTZ

We generated expression profiles of fly heads at four time points – 0 hr; 12 hrs of PTZ, the beginning phase; 2^nd day of PTZ, the latent phase; and 7^th day of PTZ, the behaviorally expressive phase (Figure 4). The full microarray data set has been deposited in the Gene Expression Omnibus http://www.ncbi.nlm.nih.gov/geo/ under accession series GSE7037, GSE7120 and GSE7200. Total number of analyzable genes in SAM showed uniformity across time points – 4369, 4637, 5259 and 4297, in that order. At 0 hr, no gene was detected as differentially expressed below 96% FDR. At 12 hrs, 2^nd day and 7^th day time-points, 23, 2439 and 265 genes were found to be downregulated at 22.76%, 13.74% and 23.03% FDR, in that order (for gene lists, see additional file 3). No upregulated gene was detected at above FDRs. Given small number of genes at 12 hrs, we considered these FDRs as the best compromise between uniformity across time points and acceptability in terms of incorporating false positives. We therefore used the above gene lists for all subsequent analysis. Detection of genes which are known to express at a relatively low level in fly brain, diss (CG4211), dnc (CG32498), dco (CG2048), per (CG2647), for example [46], as differentially expressed genes in PTZ time series demonstrated that our expression profiling was efficient. Further, a few genes which were used to validate microarray results by real time PCR corroborated the transcriptomic profiles (for details, see additional file 4). Notably, significant overlap of genes between adjacent time points [hypergeometric distribution; p = 5.3 × 10^-6 between 12 hrs and 2^nd day, p = 1.9 × 10^-22 between 2^nd and 7^th day] demonstrated that our expression profiling was efficient enough to capture the spectrum of transcriptomic alterations underlying chronic PTZ phase.

We first examined if 'seizure' genes already described in Drosophila show enrichment in our PTZ induced downregulated genes. Evidence suggests conserved mechanisms between Drosophila bang-sensitive mutants and human seizure disorders [20]. Molecular defects in five such mutants involve genes encoding mitochondrial ribosomal protein (CG7925), ADP/ATP translocase (CG16944), the citrate synthase (CG3861), an ethanolamine kinase (CG3525), and an aminopeptidase (CG5518). Prevalence of mitochondrial function disrupting defects in bang-sensitive mutants suggests that impaired energy metabolism may underlie the affected behavior [20]. Importantly, human seizure disorders have been linked to mutations in genes encoding pyruvate carboxylase, and pyruvate dehydrogenase, the two enzymes upstream of citrate synthase, as well as to mutations affecting various steps in the TCA cycle [20]. Our microarrays had all the above bang-sensitive genes spotted except CG7925. Notably, all four spotted seizure genes were downregulated by PTZ. This overlap was statistically significant (hypergeometric distribution, p = 0.004).

For providing direction to the transcriptomic changes, we first overlaid all 2574 genes present in 12 hrs, 2^nd day and 7^th day sets onto protein interactome map of D. melanogaster (interaction database BioGRID v 2.0, http://www.thebiogrid.org/index.php visualization software platform Osprey v. 1.2.0, http://biodata.mshri.on.ca/osprey/servlet/Index) which consisted of 7358 genes (vertices) and 24984 connections (edges) [47]. Out of 2574 PTZ regulated genes, 1633 were present in the fly interactome. We retrieved PTZ gene specific subinteractome that consisted of 1633 vertices and 1485 edges. In this subinteractome map, 632 genes were found to be loner, i.e., they did not show any within group interaction. We hereafter focused on within group interaction network which was comprised of 1001 genes and 1485 interactions (for interacting genes, see additional file 5). We next examined which of the 12 hrs gene(s) shows statistically significant overinteraction in the network of 1001 vertices and 1485 edges. Among 13 beginning phase, i.e., 12 hrs, genes present in the network, only the C-terminal Binding Protein (CtBP, CG7583) was found to overinteract (hypergeometric distribution after Bonferroni correction, p = 0.004; Table 3; Figure 5). CG7583 and its direct interacting partners, a network comprising of 17 genes, showed overrepresentation of various GO processes which were mainly related to transcriptional regulation and cell fate determination (Table 4). Transcriptional alteration was thus identified as the most proximal effect of PTZ. To complete adding direction, we retrieved enriched processes in genes which showed interaction with interacting partners of CG7583 in the self-interacting network of our PTZ genes, and so on, till completion; of 1001 genes, 898 were consumed by this CG7583 centered network. Enrichment of cell morphogenesis and protein import into nucleus; neuron differentiation, neurogenesis, axonogenesis, axon guidance, small GTPase mediated signal transduction, actin filament organization, cell growth, transmembrane receptor protein tyrosine kinase signaling pathway, eye photoreceptor cell differentiation, regulation of cell shape and regulation of cell size; phospholipase C activation and positive regulation of hydrolase activity was observed, in that order (for complete lists of enriched GO processes, see additional file 6). Broadly, process enrichment analysis guided by protein interactome suggested that chronic PTZ sequentially downregulates transcription, neurogenesis/axonogenesis/axon guidance and phospholipase C activation. This was consistent with original time series analysis presented above (Tables 1 and 2), except that glutamate metabolism related processes enriched on 7^th day (Table 2) were not recovered by interactome guided analysis. Seventh day genes under these categories – CG7145, CG5320 (Gdh), CG14994 (Gad1), CG1743 (Gs2) and CG8430 (Got1) – are poorly connected in the entire fly interactome (range, 0–3 interactions). It is likely that non-recovery of glutamate metabolism related processes reflect poor interaction of the genes involved. Cumulatively, the above results suggested that sequential downregulation of gene expression, neurogenesis/axonogenesis/axon guidance, phospholipase C activation and glutamate metabolism characterize our fly behavioral model.

We next examined enrichment of biochemical and signaling pathways in all 338 genes which were part of enriched biological processes in CG7583 centered network, listed in additional file 6. Although no pathway was found to be enriched after Bonferroni correction for multiple hypotheses testing, our gene set overrepresented various pathways based on unadjusted p values (Table 5). Since enrichment of these pathways is consistent with process enrichment analysis presented above – TGF-beta signaling pathway and Wnt signaling pathway are consistent with neurogenesis/axonogenesis/axon guidance, for example – we considered the enriched pathways as significant. Genes showing downregulation in PTZ fly model are indicated in the enriched pathways shown in Figures 6, 7, 8, 9, 10, 11. Remarkably, Wnt signaling which is the most significant pathway in our analysis has previously been implicated in epilepsy [48]. As depicted in Figures 6, Wnt signaling is directly connected to TGF-beta (Figure 7) and MAPK signaling (Figure 8) pathways. Further, one or more of these pathways are also linked to JAK-STAT (Figure 9) and Cell Communication (Figure 10) pathways. Therefore, except Dorso-Ventral Axis Formation (Figure 11) which does not show any direct or indirect connection with Wnt signaling, all other five pathways enriched in our gene set are marked by biological plausibility in epileptogenesis.

We further examined the relevance of fly model in epilepsy and related neurological and psychiatric conditions through mining reported transcriptomic and proteomic data related to available animal models and human patients. For this, we retrieved rat, mouse and human homologs (FLIGHT, homologene; http://www.flight.licr.org/) of above mentioned 1001 genes which were downregulated by PTZ and were part of self-interacting network (for homolog gene lists, see additional file 7). The homologs were then examined to find out if they overrepresent genes/proteins reportedly downregulated in transcriptomic and proteomic studies on rodent models or human epilepsy (for gene/protein lists, see additional file 8). Comparison with individual studies did not reveal statistically significant overlap. Since overlap was also not observed even among reported studies, we pooled down- and up-regulated genes therein separately and matched with our fly gene list. Notably, downregulated genes in fly showed significant match with reported downregulated (p = 0.017) but not upregulated (p = 0.082) genes. This direction specific overlap demonstrated that our fly model is relevant in understanding epileptogenesis.