Transcriptomic analysis in a Drosophila model identifies previously implicated and novel pathways in the therapeutic mechanism in neuropsychiatric disorders.

We have taken advantage of a newly described Drosophila model to gain insights into the potential mechanism of antiepileptic drugs (AEDs), a group of drugs that are widely used in the treatment of several neurological and psychiatric conditions besides epilepsy. In the recently described Drosophila model that is inspired by pentylenetetrazole (PTZ) induced kindling epileptogenesis in rodents, chronic PTZ treatment for 7 days causes a decreased climbing speed and an altered CNS transcriptome, with the latter mimicking gene expression alterations reported in epileptogenesis. In the model, an increased climbing speed is further observed 7 days after withdrawal from chronic PTZ. We used this post-PTZ withdrawal regime to identify potential AED mechanism. In this regime, treatment with each of the five AEDs tested, namely, ethosuximide, gabapentin, vigabatrin, sodium valproate, and levetiracetam, resulted in rescuing of the altered climbing behavior. The AEDs also normalized PTZ withdrawal induced transcriptomic perturbation in fly heads; whereas AED untreated flies showed a large number of up- and down-regulated genes which were enriched in several processes including gene expression and cell communication, the AED treated flies showed differential expression of only a small number of genes that did not enrich gene expression and cell communication processes. Gene expression and cell communication related upregulated genes in AED untreated flies overrepresented several pathways - spliceosome, RNA degradation, and ribosome in the former category, and inositol phosphate metabolism, phosphatidylinositol signaling, endocytosis, and hedgehog signaling in the latter. Transcriptome remodeling effect of AEDs was overall confirmed by microarray clustering that clearly separated the profiles of AED treated and untreated flies. Besides being consistent with previously implicated pathways, our results provide evidence for a role of other pathways in psychiatric drug mechanism. Overall, we provide an amenable model to understand neuropsychiatric mechanism in cellular and molecular terms.

Introduction

Drugs used in the treatment of psychiatric disorders are mostly known to target neurotransmitter receptors. This receptor mechanism alone however does not provide simple mechanistic interpretations for their long term clinical efficacy (Molteni et al., 2009; Zhou et al., 2009). In addition to receptor mediated acute biochemical effects that may explain short-term clinical response, these drugs are considered to exert other long term therapeutic effects that may not be directly related to receptor mechanisms (Molteni et al., 2009). These long term neuroprotective mechanisms underlying psychiatric drug action are however poorly understood (McLoughlin et al., 2009). It has been suggested that drug induced changes in gene and protein expression may ultimately translate into the overall neuroprotection. Given this, genome level expression analysis is considered to offer a promising approach to identify genes and pathways underlying neuropsychiatric conditions and mechanisms of drug action (Altar et al., 2009). For example, a meta-analysis has recently revealed that gene expression profiles of brains from persons with major depressive disorder show decreased expression of genes related to glutamate transport and metabolism, neurotrophic signaling and MAP kinase pathways, and that genes in these pathways show increased expression in the brains of rodents exposed to antidepressant treatments (Altar et al., 2009).

As with other psychiatric drugs, the long term mechanism of antiepileptic drugs (AEDs), which are also used in the treatment of various psychiatric conditions besides epilepsy, is poorly understood (Rogawski and Loscher, 2004; Johannessen, 2008; Kuzniecky et al., 2008; Nagarkatti et al., 2008; Nalivaeva et al., 2009). Conceptually, transcriptomic analysis in established mammalian models of epilepsy and AED testing can be used to gain insights into the mechanisms of action of these drugs. However, inherent complexity of mammalian brain does not render these established models as amenable to systems modeling (Gorter et al., 2006). Under these circumstances, the genetically tractable model organism Drosophila, because of its amenability to a wide variety of experimental approaches, including functional genomics (Chintapalli et al., 2007), may offer an attractive system to unravel AED mechanism.

An established model of epileptogenesis and AED testing, kindling in rodents involves long term brain plasticity in which recurrent activation of neural pathways through chemical or electrical means results in an increased susceptibility to evoked seizures that ultimately progresses to spontaneous seizures (Goddard et al., 1969; Walker et al., 2002; Husum et al., 2004; Garriga-Canut et al., 2006; McNamara et al., 2006; De Smedt et al., 2007). Like epileptogenesis, kindling is known to be associated with several non-epileptic conditions such as schizophrenia-like behaviors, hyperlocomotor activity, anxiogenic response, hyperalgesia, amnesia, spatial learning and memory, and neurodegeneration (Mortazavi et al., 2005; Pavlova et al., 2006; Szyndler et al., 2006; Akula et al., 2007; Howland et al., 2007; Omrani et al., 2007; Ma and Leung, 2010). Notably, inspired by pentylenetetrazole (PTZ) kindling in rodents, we recently described a Drosophila model of chronic PTZ induced alteration in locomotor activity (Mohammad et al., 2009), a behavior that is considered relevant in understanding neuropsychiatric conditions (Yamamoto et al., 2008; Iliadi, 2009). In this fly model, 7 days of PTZ treatment and 7 days of subsequent PTZ discontinuation progressively result in a decreased and an increased speed of startle-induced climbing in Drosophila adult, in that order. The chronic PTZ regime is responsive to AEDs; flies treated with PTZ and AED combined do not exhibit altered locomotor behavior. Importantly, downregulation of genes enriched in several processes such as transcription, cell differentiation, cell communication, neurogenesis, axonogenesis, axon guidance, and glutamate metabolism, etc., characterize the fly head transcriptome in the chronic PTZ regime. Moreover, gene expression alteration in the fly model has been found to mimic that reported in established mammalian models of epileptogenesis and human epileptic patients (Mohammad et al., 2009). These findings have suggested that the fly model may potentially be used in understanding mechanisms of action of AEDs at transcriptomic level.

Here, we describe use of the fly model to gain insights into the possible mechanism of AED action. Unlike chronic PTZ regime characterized earlier at transcriptomic level (Mohammad et al., 2009), the uncharacterized post-PTZ withdrawal part of the Drosophila model offered a simpler system for identifying potential mechanism of AEDs at transcriptomic level because it is not complicated by continued presence of the GABAA antagonist. We thus used the post-PTZ regime in the present analysis. Specifically, we tested the concept that if treatment with AEDs after PTZ withdrawal is found to rescue flies from developing increased climbing speed then analyzing CNS expression profiles of flies treated and untreated with AEDs may enable identification of transcriptomic correlates of AED action.

Materials and Methods

Behavioral pharmacology

The wild type Oregon-R strain of Drosophila melanogaster was used in the analysis. Cultures were routinely maintained at 24 ± 1°C, 60% RH, and 12 h light (9 AM to 9 PM) and 12 h dark cycle, in normal food (NF) consisting of agar–agar, maize powder, brown sugar, dried yeast, and nipagin. Standard fly handling and manipulation methods were followed. Stringency required in behavioral studies was strictly adhered to at several levels including housing conditions, exposure to anesthetic agent, light intensity, etc. Three- to 4-day-old unmated male flies were used to begin treatment at 0 day time-point (Figure 1A). Final concentration of PTZ, ethosuximide (ETH), gabapentin (GBP), vigabatrin (VGB), sodium valproate (NaVP; all from Sigma-Aldrich), and levetiracetam (LEV; Levesam 500, Nicholas Piramal) in the fly medium was 8, 3.48, 16, 24, 0.33, and 5 mg/ml, in that order. Climbing speed was measured using a semi-manual method (Mohammad et al., 2009). In this method, individual flies were first familiarized in a vertically placed glass column for 90 s and then startle-induced climbing activity was recorded using a “dot/comma” method. In “dot/comma” recording, the locomotor activity of a fly was recorded by keep pressing the dot key or the comma key of a personal computer, to record a climbing or a resting fly, in that order. Using the cursor speed, the dots and commas were accordingly transformed in the activity and rest period. Climbing speed was calculated using the following formula, s = h/t, where s = climbing speed, h = height climbed in centimeter, and t = activity period in second.

Figure 1

Study design for behavioral and transcriptomic analyses. Twelve different treatment conditions for behavioral analysis are depicted (A). Climbing speed was measured at 14th day time-point. Six treatment combinations for gene expression comparison are shown (B). Expression profiles were generated at 10th and 14th day time-points. NF, normal food; ETH, ethosuximide; GBP, gabapentin; VGB, vigabatrin; NaVP, sodium valproate, LEV, levetiracetam. See text for details.

Microarray analysis

Total cellular RNA was isolated from fly heads belonging to four biological replicates. 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 Cy3 and Cy5 (Amersham) labeled antisense RNA. The Cy3 and Cy5 labeled aRNAs (control and treated) were pooled together and precipitated, washed, air-dried, and dissolved in 18 MΩ RNAase free water. A total of 48 microarrays were hybridized, four each for 10th and 14th day flies not treated with any AED and 10th and 14th day flies treated with each of the five AEDs separately (Figure 1B). The arrays used in the experiment (Canadian Drosophila Microarray Centre, Toronto) represent over 10000 unique D. melanogaster genes and are available for distribution to academic labs. Each microarray compared RNA abundance in drug exposed flies versus flies never exposed to any drug at any time, i.e., maintained throughout in NF. Out of four slides representing four biological replicates, two were dye-swaps. Each biological replicate represented RNA isolated from heads of 120 control or treated flies. These flies were collected from four vials, with 30 flies housed in each. The control and drug fed flies were always treated in parallel. The four biological replicates represented control and treated flies collected on four different days. Microarray hybridization was set-up on eight different days, four each for 10th and 14th day time-points. A single biological replicate of each of the six comparisons to be carried out for a given time-point were processed together and used for hybridization in parallel, on a single day (Figure 2). Hybridization solution contained hybridization buffer (DIG Easy Hyb, Roche), 10 mg/ml salmon testis DNA (0.05 mg/ml final concentration, Sigma), 10 mg/ml yeast tRNA (0.05 mg/ml final concentration, Sigma), and the Cy3 and Cy5 labeled product. The mixture was denatured at 65°C and applied onto microarray slides. The slides were covered by a 24 mm × 60 mm coverslip (ESCO, Portsmouth, USA). Hybridization was carried out in a hybridization chamber (Corning) at 37°C for 16 h. After hybridization, slides were submerged in a solution containing 1× SSC and 0.1% SDS at 50°C, to remove the coverslips. Slides were washed in 1× SSC and 0.1% SDS at 50°C (three times for 15 min each) and then in 1× SSC at room temperature (twice for 15 min each). Slides were further washed in 0.1× SSC for 15 min and the liquid remaining on the slides after washing was quickly removed by spinning the slides at 600 rpm for 5 min.

Figure 2

Batch structure of microarrays. Microarrays were run in eight batches. Each batch comprised of one of the four biological replicates belonging to 12 comparisons depicted in Figure 1. Fly treatment, RNA isolation, labeling, hybridization, and scanning were carried out separately for each of the eight batches.

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). 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 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 significant analysis of microarrays (SAM 3.0, Excel Add-In), under the conditions of one class response and 100 permutations (Tusher et al., 2001). Normalized log2 ratio (635/532) of four biological replicates with balanced dye-swaps was used for microarray clustering using Acuity 4.0 (Molecular Devices). Details of RNA extraction and microarray analysis have been described previously (Mohammad et al., 2009). The full microarray data set has been deposited in the Gene Expression Omnibus under accession series GSE7156, GSE10984, GSE10985, GSE10986, GSE10987, and GSE10988.

Bioinformatics

FLIGHT was used for retrieving gene symbols and IDs. The GOTool Box (Martin et al., 2004) was used to retrieve overrepresented gene ontology (GO) biological processes in up- or down-regulated genes. The GOTool Box was used under the settings, ontology, biological process; mode, all terms; reference, genome; evidence, all–all evidence; species, D. melanogaster; GO-stats. DAVID was used for examining enrichment of pathways in Kyoto encyclopedia of genes and genomes (KEGG) database using modified Fisher exact test (Dennis et al., 2003; Huang et al., 2009). Genes showing pathway enrichment in DAVID analysis were depicted in the KEGG maps for D. melanogaster.

Results

Behavioral effect of PTZ is normalized by AEDs

We first examined the behavioral pharmacology of AEDs in the fly model. In this analysis, flies were treated with PTZ for 7 days, with an AED for 3 days, and with NF for 4 days, in that order, before climbing speed was measured (Figure 1A). In parallel, flies were treated either with NF for the entire 14 days, or with PTZ for first 7 days and then by NF for rest of the period, or with NF, AED, and NF for 7, 3, and 4 days in sequence (Figure 1A). Previously, a control climbing assay using the same method followed here detected no significant variation among various batches of NF flies housed in different vials, neither in one-way ANOVA nor in two-tailed, pair-wise Student's t-test with p-value unadjusted for multiple testing (Mohammad et al., 2009). This demonstrated that climbing speed measurement was a robust assay and pair-wise comparison at nominal p-value is reliable. In the present analysis, we used Student's t-test with nominal p-value for pair-wise comparisons. As expected from the previous report (Mohammad et al., 2009), PTZ withdrawal alone caused an increased climbing speed compared to NF (Figure 3). However, climbing speed of flies treated with either of the AEDs except ETH did not differ significantly from that of NF control. In ETH group also, though the flies showed higher speed than NF (p = 0.013), the difference was far less significant compared to that observed between PTZ withdrawal alone and NF (p = 0.00000017). In AED alone group, none except VGB caused a significantly altered speed. The observed rescuing effect of GBP, NaVP, and LEV in post-PTZ withdrawal flies was thus found not to be confounded by their locomotor effect in normally grown flies. Cumulatively, the AEDs in general normalized the behavioral perturbation caused by PTZ withdrawal.

Figure 3

Behavioral pharmacology of AEDs. Mean ± SE (n = 24) of climbing speed of flies treated and untreated with AED after PTZ withdrawal, and treated with AED alone. Note that climbing speed of flies treated with an AED after PTZ withdrawal is either insignificantly or less significantly different from the control (NF) group, compared to flies not treated with an AED after PTZ withdrawal. PTZ indicates no AED treatment after PTZ withdrawal, AED abbreviation indicates AED alone treatment, PTZ + AED abbreviation indicates AED treatment after PTZ withdrawal. Speed in the control (NF) group was compared with various treatment groups, in pair-wise comparisons. *Indicates nominal p-value. See text for details.

Transcriptomic effect of PTZ is normalized by AEDs

Differentially expressed genes

We next asked the question if AEDs in general normalize the transcriptomic perturbation caused by PTZ withdrawal. To explore this, microarray gene expression profiles of fly heads were generated at two time-points – 3rd and 7th day after PTZ withdrawal, i.e., on 10th and 14th day from the start of PTZ treatment. Twelve sets of flies were profiled, six each for 10th and 14th day time-points (Figure 2). In one set each for the two time-points, flies were not treated with any AED after PTZ withdrawal (Figure 2). In five other sets for each time-points, flies were treated with PTZ for 7 days, with one of the five AEDs for 3 days (10th day time-point), or with NF for 4 days (14th day time-point), in that order (Figure 2). Flies treated in parallel with NF throughout were used for comparison, in each of the 12 sets of microarrays. Four biological replicates comprised each set of microarray; one set compared PTZ withdrawal with NF control, and the rest compared ETH, GBP, VGB, NaVP, or LEV treatment after PTZ withdrawal with NF control (Figure 2). In a preliminary analysis, we observed an increasing enrichment of GO biological processes in differentially expressed genes up to 15% false discovery rate (FDR). Previously, a control microarray experiment that used the same method which was followed here compared NF versus NF flies and reported no differentially expressed gene below 96% FDR (Mohammad et al., 2009). Considering the above, we used 15% FDR cut-off for identifying differentially expressed genes. Genes were found to be differentially expressed in all the six sets of microarrays, in both 10th and 14th day time-points. The up- and down-regulated genes are listed in Table S1 in Supplementary Material. The numbers of differentially expressed genes and the total analyzable spots in significant analysis of microarrays (SAM) are provided in Table 1. Differentially expressed genes in AED untreated flies (Table S1 in Supplementary Material) showed significant overlap between 10th and 14th day, in a direction-specific manner: 20 genes were common between 929 upregulated genes on 10th day and 48 upregulated genes on 14th day; 10 genes were common between 49 downregulated genes on 10th day and 158 downregulated genes on 14th day. Considering the total number of unique genes represented in the arrays as 10500, these overlaps were extremely significant (hypergeometric distribution p = 9E-10 and 2E-09, in that order). These results were not surprising because a significant overlap is expected at adjacent time-points. Very small number of differentially expressed genes (or no such genes) precluded similar analysis for AED treated profiles. Overall, the above analysis proved the robustness of the array data.

Table 1

Numbers of SAM analyzable spot IDs and differentially expressed genes in microarray profiles of flies treated with or without AEDs after PTZ withdrawal.

	No AED	ETH	GBP	VGB	NaVP	LEV
10TH DAY
Analyzable spots	7877	9450	8593	7443	8775	2760
Upregulated genes	929	1	0	0	0	4
Downregulated genes	49	0	8	283	1	4
14TH DAY
Analyzable spots	6353	5107	4473	5609	5505	2651
Upregulated genes	48	0	0	42	0	0
Downregulated genes	158	7	203	648	104	5

Using 2 × 2 chi-square test with Yates’ correction for continuity, we compared the number of analyzable spots and the number of up- and down-regulated genes in flies without AED treatment with corresponding numbers in AED treated flies. A pair-wise comparison revealed extremely significant difference in these numbers, for 10th as well as 14th day time-points. The p-values obtained for all the 10 comparisons, five each for 10th and 14th day time-points, were in the range of 0.00–0.0009. All comparisons except that with 14th day GBP and 14th day VGB showed reduction in the number of differentially expressed genes by AEDs. Compared to flies without AED treatment, the GBP and VGB treated flies displayed higher number of differentially expressed genes on 14th day. Together, AED treatment in general was found to reduce the number of differentially expressed genes on 10th day, i.e., the time-point till which flies were treated with the AEDs.

Process enrichment analysis

We next examined enrichment of GO biological processes in the differentially expressed genes. In AED untreated flies, the genes upregulated on 10th day showed enrichment for GO processes related to metabolism, gene expression and cell communication, and the downregulated genes that for cellular metabolic process (Table 2). Neither upregulated nor downregulated gene set representing 14th day time-point showed enrichment for any process. In AED treated flies, differentially expressed genes in 10th day showed enrichment of processes only in the flies treated with GBP and LEV, not ETH, VGB, and NaVP (Table 2). These enriched processes were related to metabolism. In 14th day time-point, genes differentially regulated by AEDs were found to enrich various processes. These processes were related to metabolism, response to toxin, skeletal myofibril assembly, spindle elongation, translation, CMP-sialic acid transport, etc. (Table 2). Notably, gene expression and several related processes such as transcription, transcription from RNA polymerase II promoter, RNA metabolic process, RNA biosynthetic process etc. were enriched in genes upregulated in flies that were not treated with any AED (Table 2). Also notable was the enrichment of cell communication, besides others, in genes upregulated in AED untreated group (Table 2). Combined with the previous observation (Table 1) that AEDs in general normalize the transcriptomic perturbation, the GO enrichment analysis further suggested that AEDs cause normalization of perturbation in several processes including gene expression and cell communication.

Table 2

Enriched GO processes in differentially expressed genes in flies treated with or without AEDs after PTZ withdrawal.

GO_ID	Term	p-value*
10TH DAY, UPREGULATED, PTZ WITHDRAWAL, NO AED
GO:0008152	Metabolic process	3.53E-09
GO:0044238	Primary metabolic process	3.96E-09
GO:0006139	Nucleobase, nucleoside, nucleotide, and nucleic acid metabolic process	2.48E-07
GO:0043283	Biopolymer metabolic process	6.22E-07
GO:0044237	Cellular metabolic process	1.28E-06
GO:0043170	Macromolecule metabolic process	1.58E-06
GO:0009987	Cellular process	1.06E-05
GO:0034960	Cellular biopolymer metabolic process	2.23E-05
GO:0009058	Biosynthetic process	5.58E-05
GO:0044249	Cellular biosynthetic process	6.06E-05
GO:0044260	Cellular macromolecule metabolic process	7.72E-05
GO:0034961	Cellular biopolymer biosynthetic process	0.0001
GO:0043284	Biopolymer biosynthetic process	0.0001
GO:0006350	Transcription	0.00016
GO:0010467	Gene expression	0.00017
GO:0065007	Biological regulation	0.0002
GO:0050794	Regulation of cellular process	0.0004
GO:0006366	Transcription from RNA polymerase II promoter	0.0004
GO:0034645	Cellular macromolecule biosynthetic process	0.0006
GO:0009059	Macromolecule biosynthetic process	0.0006
GO:0016070	RNA metabolic process	0.0008
GO:0050789	Regulation of biological process	0.001
GO:0010556	Regulation of macromolecule biosynthetic process	0.003
GO:0006357	Regulation of transcription from	0.004
	RNA polymerase II promoter
GO:0065003	Macromolecular complex assembly	0.004
GO:0031326	Regulation of cellular biosynthetic process	0.004
GO:0009889	Regulation of biosynthetic process	0.004
GO:0043933	Macromolecular complex subunit organization	0.005
GO:0045449	Regulation of transcription	0.007
GO:0006351	Transcription, DNA-dependent	0.007
GO:0032774	RNA biosynthetic process	0.008
GO:0019219	Regulation of nucleobase, nucleoside, nucleotide, and nucleic acid metabolic process	0.009
GO:0010468	Regulation of gene expression	0.012
GO:0034621	Cellular macromolecular complex subunit organization	0.017
GO:0080090	Regulation of primary metabolic process	0.021
GO:0007154	Cell communication	0.031
GO:0034622	Cellular macromolecular complex assembly	0.036
GO:0060255	Regulation of macromolecule metabolic process	0.036
10TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, NO AED
GO:0044237	Cellular metabolic process	0.024
10TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, GBP
GO:0006629	Lipid metabolic process	0.016
GO:0009636	Response to toxin	0.05
GO_ID	Term	p-value*
10TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, LEV
GO:0048252	Lauric acid metabolic process	0.002
GO:0031000	Response to caffeine	0.007
GO:0014074	Response to purine	0.007
GO:0009404	Toxin metabolic process	0.02
GO:0017143	Insecticide metabolic process	0.02
GO:0006805	Xenobiotic metabolic process	0.02
GO:0009410	Response to xenobiotic stimulus	0.02
GO:0014070	Response to organic cyclic substance	0.04
GO:0043279	Response to alkaloid	0.04
GO:0017085	Response to insecticide	0.046
14TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, ETH
GO:0009636	Response to toxin	0.0051
14TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, GBP
GO:0008152	Metabolic process	3.36E-06
GO:0016052	Carbohydrate catabolic process	0.00012
GO:0046164	Alcohol catabolic process	0.0024
GO:0050896	Response to stimulus	0.0053
GO:0006091	Generation of precursor metabolites and energy	0.006
GO:0044237	Cellular metabolic process	0.006
GO:0009636	Response to toxin	0.02
GO:0044275	Cellular carbohydrate catabolic process	0.03
GO:0044248	Cellular catabolic process	0.03
14TH DAY, UPREGULATED, PTZ WITHDRAWAL, VGB
GO:0014866	Skeletal myofibril assembly	0.047
14TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, VGB
GO:0008152	Metabolic process	6.04E-11
GO:0000022	Mitotic spindle elongation	5.01E-07
GO:0051231	Spindle elongation	6.31E-07
GO:0044237	Cellular metabolic process	1.37E-06
GO:0044238	Primary metabolic process	3.30E-05
GO:0019538	Protein metabolic process	6.66E-05
GO:0044267	Cellular protein metabolic process	9.63E-05
GO:0044248	Cellular catabolic process	0.0006
GO:0006412	translation	0.0007
GO:0007052	Mitotic spindle organization	0.001
GO:0050896	Response to stimulus	0.002
GO:0006950	Response to stress	0.005
GO:0009056	Catabolic process	0.009
GO:0007051	Spindle organization	0.03
14TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, NAVP
GO:0009636	Response to toxin	2.86E-06
GO:0050896	Response to stimulus	0.0003
GO:0009056	Catabolic process	0.0003
GO:0016052	Carbohydrate catabolic process	0.004
GO:0005975	Carbohydrate metabolic process	0.008
GO:0008152	Metabolic process	0.03
14TH DAY, DOWNREGULATED, PTZ WITHDRAWAL, LEV
GO:0015782	CMP-sialic acid transport	0.03
GO:0015789	UDP-N-acetylgalactosamine transport	0.03
GO:0015757	Galactose transport	0.03

Gene expression and cell communication related processes are highlighted in bold.

*After Bonferroni correction.

Pathway enrichment analysis

To gain further insights into the normalizing effect of AEDs on PTZ withdrawal induced transcriptomic perturbation, we next examined if gene expression (GO:0010467) and cell communication (GO:0007154) related upregulated genes in AED untreated flies on 10th day enrich specific pathways. Notably, several pathways were found to be overrepresented in the two sets of genes (Table 3) – the gene set belonging to gene expression category showed enrichment for spliceosome, RNA degradation and ribosome, and the gene set belonging to cell communication enriched inositol phosphate metabolism, phosphatidylinositol signaling, endocytosis, and hedgehog signaling. The upregulated genes are mapped on to KEGG pathways for visualization (Figures 4–10).

Table 3

Enriched processes in gene expression and cell communication related upregulated genes in AED untreated flies.

Term	p-value*
GENE EXPRESSION
dme03040:spliceosome	0.001
dme03018:RNA degradation	0.006
dme03010:ribosome	0.049
CELL COMMUNICATION
dme00562:inositol phosphate metabolism	0.012
dme04070:phosphatidylinositol signaling system	0.02
dme04144:endocytosis	0.027
dme04340:hedgehog signaling pathway	0.046

*Nominal.

Figure 4

Genes mapping to KEGG pathway for spliceosome. Pink boxes represent upregulated genes in AED untreated flies. These genes are Hsp68 (Heat shock protein 68), U2af50 (U2 small nuclear riboprotein auxiliary factor 50), snRNP69D (small nuclear ribonucleoprotein at 69D), and snRNP70K (small nuclear ribonucleoprotein 70K). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Figure 10

Genes mapping to KEGG pathway for hedgehog signaling. Pink boxes represent upregulated genes in AED untreated flies. These genes are Rab23 (CG2108 gene product from transcript CG2108-RA), dpp (decapentaplegic), and fu (fused). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Genes mapping to KEGG pathway for RNA degradation. Pink boxes represent upregulated genes in AED untreated flies. These genes are Csl4 (CG6249 gene product from transcript CG6249-RA), Dis3 (CG6413 gene product from transcript CG6413-RA), and Rrp6 (CG7292 gene product from transcript CG7292-RB). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Genes mapping to KEGG pathway for ribosome. Pink boxes represent upregulated genes in AED untreated flies. These genes are RpL15 (Ribosomal protein L15), RpL28 (Ribosomal protein L28), RpL38 (Ribosomal protein L38), RpL39 (Ribosomal protein L39), and RpL5 (Ribosomal protein L5). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Genes mapping to KEGG pathway for inositol phosphate metabolism. Pink boxes represent upregulated genes in AED untreated flies. These genes are IP3K2 (Inositol 1,4,5-triphosphate kinase 2), Ipk1 (CG30295 gene product from transcript CG30295-RB), Ipp (Inositol polyphosphate 1-phosphatase), and PIP5K59B (CG3682 gene product from transcript CG3682-RA). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Genes mapping to KEGG pathway for phosphatidylinositol signaling system. Pink boxes represent upregulated genes in AED untreated flies. These genes are IP3K2 (Inositol 1,4,5-triphosphate kinase 2), Ipk1 (CG30295 gene product from transcript CG30295-RB), Ipp (Inositol polyphosphate 1-phosphatase), and PIP5K59B (CG3682 gene product from transcript CG3682-RA). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Genes mapping to KEGG pathway for endocytosis. Pink boxes represent upregulated genes in AED untreated flies. These genes are Gap69C (GTPase-activating protein 69C), Hsp68 (Heat shock protein 68), PIP5K59B (CG3682 gene product from transcript CG3682-RA), TSG101 (tumor suppressor protein 101), Vps20 (Vacuolar protein sorting 20), alpha-Adaptin (CG4260 gene product from transcript CG4260-RA), cenG1A (centaurin gamma 1A), and schizo (K12495 IQ motif and SEC7 domain-containing protein). Green boxes represent other members in the pathway database for Drosophila melanogaster.

Microarray clustering

Transcriptomic analysis so far depended exclusively on differentially expressed genes selected using a FDR cut-off. To further examine if AEDs indeed remodel PTZ withdrawal induced transcriptomic perturbation, we next clustered the microarray profiles of flies with or without AED treatment. The hierarchical clustering that uses expression values of all the analyzable spots in the microarrays clearly separated the time series profiles of AED treated and untreated flies (Figure 11). Microarray clustering thus overall confirmed the transcriptome remodeling effect of AEDs.

Figure 11

Hierarchical clustering of microarrays. Microarrays of 10th and 14th day time-points are clustered to examine separation of profiles of flies untreated with an AED and flies treated with the AEDs ETH (A), NaVP (B), LEV (C), VGB (D), and GBP (E) after PTZ withdrawal. Note clear separation of the two groups. The AED untreated group is indicated by PTZ, and the AED treated groups by AED abbreviation. Jaccard similarity metric and average linkage methods were used for clustering.

Discussion

We have shown here that AEDs normalize long term behavioral and transcriptomic alterations induced by PTZ withdrawal in a Drosophila model. Our evidence further suggests that AEDs’ transcriptomic effect is mediated by neutralization of upregulated genes related to various processes including gene expression and cell communication. Furthermore, our results suggest that AEDs neutralize upregulation of genes belonging to several pathways. In the cell communication category, these pathways are inositol phosphate metabolism, phosphatidylinositol signaling, endocytosis, and hedgehog signaling. In the gene expression category, these pathways are spliceosome, RNA degradation, and ribosome. It is important to note here that cell communication pathways have previously been implicated in the therapeutic mechanisms of AEDs and antipsychotic drugs in diverse studies. For example, biochemical and neurobiological evidence has earlier suggested that NaVP and LEV inhibit inositol metabolism and/or phosphatidylinositol signaling (Simister et al., 2007; Xu et al., 2007; Nagarkatti et al., 2008; Tokuoka et al., 2008; Teo et al., 2009; Yamamura et al., 2009). Also, cell biological evidence has shown NaVP induced reduction in endocytosis, a process that is linked to phosphatidylinositol signaling (Xu et al., 2007). Similarly, transcriptomic evidence showing increased expression of endocytosis related genes in phenytoin resistant kindled rats, a model of epileptogenesis, has previously been presented (Zeng et al., 2009). Further, several antipsychotic drugs have recently been found to regulate hedgehog signaling (Lauth et al., 2010), a pathway that, besides its established role in adult CNS, is also known for its growth enhancing effect in the adult brain (Tsuboi and Shults, 2002; Dellovade et al., 2006).

Regarding gene expression related pathways, though transcription factors and mRNA expression have earlier been implicated (Atmaca, 2009; Heinrich et al., 2009; Christensen et al., 2010; Girgenti et al., 2010), the involvement of spliceosome, RNA degradation, and ribosomal pathways has not been well documented in the therapeutic mechanism of AEDs. It is therefore notable that our study has identified these three pathways as additional candidates in the therapeutic mechanism. Altered regulation of RNA metabolism including splicing, mRNA stability, etc., are known to be associated with various neurological and psychiatric disorders (Licatalosi and Darnell, 2006; Anthony and Gallo, 2010). Also, recent biochemical studies have shown that phosphatidylinositol signaling is linked to mRNA processing and translation (Kiefer et al., 2009; Laserna et al., 2009; Lewis et al., 2011). The above evidence supports our results that implicate a spectrum of gene expression and cell communication pathways in the therapeutic mechanism of neuropsychiatric disorders. It is intriguing though that we did not find enrichment of processes related to synaptic function in differentially expressed genes, as would be expected for proconvulsant and anticonvulsant drugs. It may however be noted here that our transcriptomic results are based on the analysis of only a subset of the total genes represented on the microarrays, due to the loss of a large number of genes during downstream processing of the data. More sensitive labeling, hybridization, and detection system, such as those provided by Affymetrix platform, could have possibly minimized this loss and resulted in the identification of additional differentially expressed genes and, in turn, enrichment of additional processes and pathways.

Acute biochemical effects of drugs used in treating neuropsychiatric disorders do not provide simple mechanistic interpretations for the observed neuroprotection (Zhou et al., 2009). Instead, genomic effects of these drugs are considered to ultimately translate into the overall neuroprotection (Altar et al., 2009). Given this, it is significant that transcriptomic analysis in a novel Drosophila model has provided evidence that is consistent with known effect of drugs reported in mammalian studies. Amenable to various experimental approaches, the Drosophila model may thus offer a unique opportunity to further understand the psychotropic drug mechanism in cellular and molecular terms. Methods in fly genetics, for example, may be used to functionally validate candidate genes in the pathways. Similarly, epigenetic approaches could be applied to understand long term drug effects in brain plasticity. Besides, the fly model may directly be used for drug screening using behavioral and functional genomic readouts. It is tempting to note here that we are currently implementing the fly model in screening drugs approved by Food and Drug Administration (USA) toward potential repurposing. The preliminary results obtained so far are indeed encouraging. Further screening of positive compounds in a battery of rodent models available for epileptogenesis, and neuropsychiatric and neurodegenerative disorders will be crucial. In brief, the fly model promises to accelerate the pace of discovery in the area of CNS disorders and therapy.

Supplementary Material

The Supplementary Material for this article can be found online at http://frontiersin.org/neuroscience/neurogenomics/paper/10.3389/fnins.2011.00161/abstract

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.