284 Quantitative RT-PCR mRNA was reverse transcribed with random hexamers to generate single-stranded cDNA transcripts according to the manufacturer’s protocol of the High Capacity cDNA Ar- chive Kit (Applied Biosystems). For qPCR reactions, cDNA from three individuals was taken for each group. A prede- signed TaqMan Low Density Array (TLDA) in a 96-well format (Applied Biosystems) was used, and duplicates were performed for the expression of each gene in each sample (for the accession numbers of the probes used, see Supple- mentary file 1). A total amount of 100 ng of cDNA for each sample, RNase-free water, and 50 μl of TaqMan Universal PCR Master Mix were loaded into each reservoir to a total volume of 100 μl/reservoir. Real-time PCR amplifications were performed using a 7900 HT Real-Time PCR Sys- tem (Applied Biosystems). The thermocycling condi- tions for the reactions were a 10-min initial setup at 95°C, followed by 40 cycles of 30 s denaturation at 97°C and 1 min annealing/extension at 60°C. The results for each plate were analyzed using ABI PRISM SDS2.3 (Applied Biosystems) software to calculate the raw value of Cq (quantification cycle, also known as the threshold cycle, crossing point, or take-off point) in each well. 18 S rRNA was used as an endogenous control, and Rpl32 and Ppia were used as reference genes. The levels of expressed genes were measured using a relative quantitative method (Pfaffl 2001), with the Ppia and Rpl32 genes as reference genes and a control group (saline treated) as a calibrator (Applied Biosystems, UserBulletin#2). Data analysis Genes whose Cq values were higher than 38 or undetermined in more than a half of the samples per analyzed brain region were excluded from analysis (for each structure and each drug, there were 10 time groups, n06 per group). Statistical analysis was performed on the list of probes detected in each structure after each specific drug treatment. All statistical analyses of the qPCR data used ΔCq values (Cq normalized to endoge- nous controls), whereas graphs represented the fold change (2 −ΔΔCq) ± SEM. Data were analyzed by two-way analysis of variance (ANOVA) with the treatment type-specific factor (single dose or repeated) and time factor. Correction for mul- tiple testing was applied separately for each of the main factors from the ANOVA by controlling for the percentage FDR (Benjamini and Hochberg 1995). ANOVA was followed by a Bonferroni test as a post hoc analysis to determine where significant differences occurred. The threshold level for sta- tistical significance was p<0.05. Psychopharmacology (2012) 223:281–297 Gene ontology analysis The functional annotation analysis tool DAVID 2008 was used to identify over-represented gene ontology (GO) groups among the gene expression profiles and to group the genes into functional classifications (Huang et al. 2009). The lists of gene symbols (Entrez IDs), together with a particular background annotation distribution (sets of genes for which mRNA levels were detected in each brain struc- ture), were used as the input files. Over-represented GO t e r m s ( G O T E R M _ C C _ A L L , G O T E R M _ B P _ A L L , GOTERM_MF_ALL) were defined as having at least three transcripts and a maximum EASE score of 0.1. Identification of transcription factor binding sites enriched in coregulated transcripts Identification of over-represented transcription factor binding sites (TFBSs) was performed using the cREMaG database (Piechota et al. 2010) by submitting Entrez IDs lists of genes that reached statistical significance for each ANOVA factor. Functional promoter sequences were identified by carrying out alignments between the 5′ upstream sequences of orthol- ogous rat and human genes. The analysis was performed with the default parameters of a 70 % conservation threshold and a maximum number of 10 most conserved TFBSs in noncoding regions between 5,000 bp upstream and 1,000 bp downstream of the transcriptional start site. SVD analysis Data on multiple gene expression at sequential time points were analyzed using the SVD as a means to capture dominant trends, called characteristic modes. SVD was performed for the data, which were prepared as follows: first, the whole dataset was subjected to ANOVA (Brandt 1997). The aim of this procedure was to select a set of genes for which the changes during our experiments were statistically significant. For each time series, the value Fc was calculated: s2A Fc ¼ s2 W where: s2A ¼ t and s2W ¼ n 1 1 t X niðxi 1 i¼1 xÞ2; t ni X X xij t i¼1 j¼1 x 2 i