INTRODUCTION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS WEB SITE REFERENCES REFERENCES

The molecular basis for many neuropsychiatric disorders remains obscure (Owen et al. 2000). These diseases frequently have important genetic contributions, but it has been difficult to identify the relevant genes because of the complexities of human genetic analyses. Importantly, the neuroanatomical regions responsible for the deficits of the neuropsychiatric disorders are often uncertain. Vivid insights into the pathogenesis of these disorders could be obtained if it were possible to obtain an extensive sampling of gene expression patterns in three dimensions for both normal and diseased specimens. In unicellular systems, useful understanding of gene networks have been obtained from high-throughput gene expression methodologies, exemplified by microarrays (Brown and Botstein 1999), gene chips (Lipshutz et al. 1999), and SAGE (Velculescu et al. 1995). Nevertheless, these valuable techniques have yet to be systematically applied to understanding how the three-dimensional (3D) structure of multicellular organisms is constructed by their genomes. Classical technologies, such as in situ hybridization (ISH) or immunohistochemistry, give high resolution images of gene expression within the brain, but are low-throughput procedures, making it difficult to obtain a representative survey of the genome under a variety of experimental situations. It is possible to image metazoan gene expression in vivo, but currently, these technologies only permit, at most, the examination of a few genes at a time (Gambhir et al. 1999; Louie et al. 2000; Zacharias et al. 2000).

Here, a method called voxelation is used to investigate the gene expression changes that occur in the mouse brain as a result of a pharmacological model of Parkinson's disease (PD). The name of the method is derived from the term voxel, which is a cubic 3D image element. Conceptually, voxelation is simple, and uses the direct creation of voxels (cubes) spatially registered with the brain. RNA extracted from the voxels is analyzed using high-throughput techniques, allowing 3D gene expression patterns to be deduced. Voxelation thus results in multiple volumetric maps of gene expression in the brain, similar to the images reconstructed using biomedical imaging systems, such as PET and CT. The essence of the idea is the simplification of complex 3D anatomy into arrays of biochemical samples, facilitating a high-throughput analysis.

	RESULTS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS WEB SITE REFERENCES REFERENCES

A Mouse Model of PD

The PD model was created by administration of toxic doses of methamphetamine (MA) to C57BL/6J mice. At the doses used, the MA model of PD has been reported to cause substantial loss (45%) of tyrosine hydroxylase-positive dopaminergic cells in the substantia nigra, as well as destruction of dopaminergic nerve terminals in the neostriatum (Sonsalla et al. 1996). A more widely used pharmacological model of PD, the MPTP model, results in similar neuropathological effects.

To confirm the induction of the PD phenotype by MA, the striata of treated mice and untreated controls were assayed for dopamine (DA), and its metabolites dihydroxyphenylacetic acid (DOPAC) and homovanillic acid (HVA; Fig. 1A-C). Statistically significant decreases in all three analytes were found in the striata of MA-treated mice. In addition, levels of tyrosine hydroxylase (TH), the rate-limiting enzyme for biosynthesis of dopamine, were assessed in the substantia nigra using real-time quantitative RT-PCR (QRT-PCR; Fig. 1D). Consistent with the reported loss of 45% of TH-positive dopaminergic neurons in the substantia nigra (Sonsalla et al. 1996), the real-time QRT-PCR studies found a statistically significant decrease of 66% in TH mRNA. The specificity of changes in DA and its metabolites was assessed by quantitating levels of a distinct neurotransmitter system. No significant changes were found in striatal levels of 5-hydroxytryptamine (5-HT, or serotonin) and its metabolite 5-hydroxyindoleacetic acid (5-HIAA; Fig. 1E,F).

View larger version (30K): [in this window] [in a new window]   Figure 1   Neurotransmitter and metabolite changes in the MA-treated mice. For all histograms, except D, n = 4 controls, 4 MA-treated. The histograms show mean ± SEM, (**) p < 0.01, (*) p < 0.05. (A) Changes in striatal dopamine. One-tailed t-test, t = 2.53, df = 6, p = 0.02. (B) Changes in striatal dihydroxyphenylacetic acid (DOPAC). One-tailed t-test, t = 3.97, df = 6, p = 0.004. (C) Changes in striatal homovanillic acid (HVA). One-tailed t-test, t = 3.84, df = 6, p = 0.004. (D) Changes in substantia nigra tyrosine hydroxylase (TH). One-tailed t-test, t = 6.18, df = 2, p = 0.013. The peak levels of substantia nigra TH from the normal and MA brains, left and right, were estimated from real-time QRT-PCR analysis of voxels G3, G4, H3, and H4 (Fig. 4B). n = 2 controls (left and right), 2 MA-treated (left and right). (E) Changes in striatal 5-hydroxytryptamine (5-HT). Two-tailed t-test, t = 0.139, df = 6, p = 0.89. (F) Changes in striatal 5-hydroxyindoleacetic acid (5-HIAA). Two-tailed t-test, t = 1.44, df = 6, p = 0.20.

Gene Expression Relationships

Brains from control and MA-treated mice were divided into 40 voxels by slicing each brain into 10 coronal sections, and cutting each of the slices into four voxels, consisting of superior and inferior, left and right (Fig. 2). Each of the 40 voxels was then analyzed using a 9000-gene microarray. To explore the relationships between brain regions, we first examined the correlations between each of voxels in terms of gene expression levels. A subset of the entire 9000-gene data set was used for this analysis, consisting of those genes most strongly differentially expressed between the anterior half (20 voxels) and posterior half (20 voxels) of the normal brain. To identify these genes, we chose the outliers with p < 0.05 based on a T statistic, resulting in a total of 1189 genes. The rationale here was that the large differences may indicate interesting genes involved in brain development.

View larger version (31K): [in this window] [in a new window]   Figure 2   Voxelation scheme. The slices are labeled A through J, from anterior to posterior. Slice A corresponds to section 8 (6.94 mm interaural, 3.14 mm bregma) of the Mouse Brain Library (Rosen et al. 2000; Williams 2000; http://www.nervenet.org/MBL/mbl.html), slice B to section 11 (5.78 mm interaural, 1.98 mm bregma), slice C to section 12 (5.50 mm interaural, 1.70 mm bregma), slice D to section 14 (4.39 mm interaural, 0.74 mm bregma), slice E to section 16 (3.80 mm interaural, 0.00 mm bregma), slice F to section 17 (2.98 mm interaural, 0.82 mm bregma), slice G to section 20 (1.68 mm interaural, 2.12 mm bregma), slice H to section 24 (0.00 mm interaural, 3.80 mm bregma), slice I to section 28 (1.31 mm interaural, 5.09 mm bregma), slice J to section 32 (2.44 mm interaural, 6.24 mm bregma). The anterior-posterior coordinates of the Mouse Brain Library (interaural and bregma relative distances) are as described (Franklin and Paxinos 1997). Odd-numbered voxels (e.g., A1, A3, B1) are from the right side of the brain, and even-numbered voxels (e.g., A2, A4, B2) are from the left.

The resulting voxel × voxel kinships are shown using spatial correlation matrices in Figure 3A. To assess the replicability of the voxelation strategy, the data were divided into separate spatial correlation matrices for the left and right halves of both the normal and MA brains. Because there are 20 voxels for each half of the brain, 20 × 20 cross-correlation matrices were obtained, which are shown as 20 × 20 voxel images, with the matrix elements intensity-color-coded to show the correlations (which can range from +1 to 1). Each of the relevant images includes a color bar to show correlation values. The rows and columns in the matrices are numbered from 1 to 20, which follow the order in which voxels were harvested from the brains. Hence, voxels 1, 2, 3, ... from the left half correspond to voxels A2, A4, B2, ..., whereas voxels 1, 2, 3, ... from the right half correspond to voxels A1, A3, B1, ... (see Fig. 2).

View larger version (79K): [in this window] [in a new window]   Figure 3   Correlation matrices. (A) Spatial (voxel × voxel) correlation analyses using genes significantly differentially expressed between the anterior and posterior halves of the brain. The matrices for the normal and PD brain are shown with the data for the left and right halves separated out. The correlation between any pair of voxels as judged by gene expression levels is read by looking along the relevant row and column, and finding the intersection. The color of the corresponding element gives the correlation between that pair of voxels by reference to the scales (right). The diagonals are the autocorrelations of the gene expression patterns for each gene and are (and should be) equal to 1. All other correlations must be between 1 and 1. The numbering is such that voxels A2, A4, ..., I2, I4 correspond to voxels 1, 2, ..., 19, 20, respectively for the left halves of the brain, and voxels A1, A2, ..., I1, I3 correspond to voxels 1, 2, ..., 19, 20 for the right halves. (B) Expression (gene × gene) correlation analyses using the 6000 most strongly expressed genes when averaged across the 40 voxels. The correlation over voxels for any pair of genes is read by looking along the relevant row and column, and finding the intersection. The genes in the normal omnibus matrix (left and right halves combined) were parsimoniously ordered based on a similarity metric. The same gene order was used to construct the corresponding matrix for the MA brain, and the data were then separated into matrices for the left and right halves of the brains. (C) Monte Carlo analysis comparing similarity of gene expression matrices shown in B for left normal with right normal, and left MA with right MA. The distributions for each brain show the discrepancy between randomly selected pairs of permuted matrices quantitated using the Frobenius norm of the resulting difference matrix. The vertical lines show the differences between the actually observed left and right matrices for each brain. For both the normal and MA brains, the similarity of the gene expression matrices between the left and right halves was highly significant (p < 0.0001 for the normal brain; p = 0.005 for the MA brain). (D) Correlation analyses identify conserved networks of highly correlated gene expression clusters. The gene expression correlation matrices show the subset of genes common to both specimens that are anterior/posterior differentially expressed and display a spatial expression correlation coefficient of >0.75 with at least one other gene within the same brain. The genes in the normal omnibus matrix (left and right halves combined) were parsimoniously ordered using a similarity metric, and the same gene order was used to construct the corresponding matrix for the MA brain. The data were then separated into matrices for the left and right halves of the brains. Two mutually exclusive clusters of coregulated genes emerged: cluster 1 (genes 1-23) and cluster 2 (genes 24-55).

Because the brain shows a high degree of bilateral symmetry, we would expect that the left and right correlation matrices should be comparable. Figure 3A shows that the left and right matrices for both the normal and MA brains are, in fact, strikingly similar. This similarity extends both within experimental groups (e.g., left normal and right normal) and between the groups (e.g., left normal and left MA). A Monte Carlo analysis, conceptually similar to one described below, showed that the similarity between the four correlation matrices was highly significant (p < 0.0001 in all cases), implying excellent reproducibility of the voxelation strategy. Another feature might also be expected from the correlation matrices: the closer together two voxels are (i.e., nearer the diagonal), the more correlated they might be anticipated to be in terms of gene expression, whereas the further apart (i.e., toward the lower left and upper right parts of the matrices), the less correlated. This expected regional kinship can be readily confirmed from the matrices.

The correlation analyses discussed so far are spatial (i.e., voxel × voxel) and provide information on the correlation between voxels based on gene expression levels. However, it is also important to know the correlations between different genesfor example, if one gene shows a strongly regulated expression pattern across voxels, how does this correlate (positively or negatively) with the expression patterns of other genes? To investigate this, we used the 40 voxels as 40 realizations of gene expression, and took a cross correlation with respect to genes to reveal the degree to which the expression of one gene is correlated with that of others across the voxels (gene × gene matrices). For the purposes of this analysis, we used the 6000 most strongly expressed genes when averaged across the 40 voxels of the normal brain to construct gene expression correlation matrices for both control and MA-treated mice (Fig. 3B). The genes in the normal matrix were parsimoniously ordered based on a similarity metric, and the same gene order was used to construct the corresponding matrix for the MA brain. As an independent assessment of the replicability of voxelation, the data were then separated into matrices for the left and right halves of the brain for both control and MA-treated groups.

Strikingly, the matrices for the left and right brain halves were very similar within the control and MA groups (i.e., comparing left normal with right normal, and left MA with right MA), as judged using Monte Carlo statistics (p < 0.0001 for the normal brain; p = 0.005 for the MA brain; Fig. 3C), showing excellent reproducibility of the voxelation strategy. Although not as visually obvious, this statistically significant similarity also extended between groups (i.e., normal left vs. MA left, p = 0.045; normal right vs. MA right, p = 0.007), indicating conservation of gene expression relationships between the control and MA-treated brains, despite the shifts in gene expression that occur as a result of the MA treatment.

To gain further insights into gene expression in the normal and MA brains, a subset of the data was extracted from the anterior/posterior differentially expressed genes, consisting of those genes with a spatial expression correlation coefficient of >0.75 with at least one other gene in both normal and MA brains. This procedure should identify networks of coregulated genes conserved between the two brains. The results of the analysis are presented in Figure 3D as gene expression correlation matrices (gene × gene; see also Table 1). Genes in the normal matrix were parsimoniously ordered based on the similarity metric used for Figure 3B, and the same gene order was used to construct the corresponding matrix for the MA brain. For both brains, the data were again separated into matrices for the left and right halves. Similar to the overall data for the 6000 most strongly expressed subset of genes (Fig. 3B), there was highly significant left/right correspondence (i.e., left normal and right normal, left MA and right MA) for the coregulated gene subset in Figure 3D (Monte Carlo simulation, p < 0.0001), confirming replicability of the voxelation strategy. This conservation was also present in comparisons between experimental groups (i.e., left normal and left MA, right normal and right MA), implying that the coregulated networks of genes are independently maintained in both the normal and MA brains.

Spatial Gene Expression Patterns

Interestingly, the sorted data in the correlation matrices of the coregulated subset shown in Figure 3D revealed two mutually exclusive clusters of genes. Cluster 1 (genes 1-23) was positively correlated within itself, and negatively correlated with cluster 2 (genes 24-55), and vice versa. The spatial map of gene expression variation across the voxels for the selected subset of genes in both the normal and MA brains is shown in Figure 4A. The figure shows that for both the normal and PD brains, cluster 1 is most strongly expressed in the anterior part of the brain, whereas cluster 2 is most strongly expressed in the posterior. The region in which cluster 2 is most strongly expressed corresponds to voxels 33-36 (I1-I4), and includes the cerebellum (Fig. 2), suggesting that the genes in cluster 2 may be particularly important in specifying this region of the brain. Figure 4A also shows that although the mutually dependent network of spatially coregulated gene clusters is maintained within each brain, the expression patterns are modified in the MA brain compared with the normal brain, both for cluster 1 and cluster 2. 

View larger version (46K): [in this window] [in a new window]   Figure 4   Spatial gene expression patterns for the subset of correlated genes. (A) Spatial expression patterns of the genes shown in Figure 3D for the normal and MA brains. The relative level of expression of any gene in any voxel is read by looking along the relevant row and column, finding the intersection, and referring to the scales. The voxel numbering in the columns of the matrix is such that voxels A1, A2, A3, ..., I2, I3, I4 correspond to voxels 1, 2, 3, ..., 37, 38, 39, 40, respectively. The genes are in the same order as for Figure 3D. The two clusters of genes are apparent, and although these have highly conserved patterns of expression within the normal and PD brains, these patterns are somewhat divergent between the two brains. (B) Level of TH expression as judged using real-time QRT-PCR. TH expression levels are found to peak in voxels corresponding to the olfactory bulb (A1-A4, voxels 1-4) and the substantia nigra (G3, G4, corresponding to voxels 27, 28; and H3, H4 corresponding to voxels 31, 32). However, the level of TH in the substantia nigra of the PD brain is substantially decreased compared with the normal brain (Fig. 1D). (C) The Nfl gene is present as two separate spots on the microarrays, corresponding to genes 39 and 50 in A. This provided an opportunity to assess within-array replicability, which was excellent (for the normal brain, r = 0.96, F[1,38] = 443.65, p < 0.0001; for MA brain, r = 0.90, F[1,38] = 153.80, p < 0.0001). (D) Expression pattern of the precerebellin-1 gene in the normal and PD brain. The line drawing shows the anatomy of the midlevel transverse section employed, which corresponds to section 8 of the Mouse Brain Library (interaural 5.40 mm, bregma 4.60 mm). (Olf) Olfactory lobes; (Str) striatum; (Hi) hippocampus; (Ce) cerebellum. Imaging software smoothed the expression patterns over the voxels. Gene expression levels can be deduced by reference to the pseudocolor scales.

Strikingly, when genes with a spatial expression correlation coefficient of >0.75 were extracted from the 6000 most strongly expressed genes rather than the anterior/posterior differentially expressed subset, a similar pattern of two gene clusters with anterior and posterior regional expression was uncovered (Supplementary Fig. 1; Supplementary Table 1 available at www.genome.org). This suggests that a fundamental property of gene expression in the brain is the distinction between anterior and posterior, including the cerebellum.

To confirm the anatomical registration of voxels from the normal and MA brains, RNA from the voxels was assayed for tyrosine hydroxylase (TH) abundance using real-time quantitative RT-PCR (QRT-PCR). TH is strongly expressed in the olfactory bulbs and the substantia nigra (Min et al. 1994), and was therefore expected to be expressed in voxels A1-A4 (voxels 1-4; olfactory bulbs) and the ventral voxels of slices G and H, that is, G3, G4 (voxels 27, 28), and H3, H4 (voxels 31, 32; substantia nigra). The real time QRT-PCR was performed using the same RNA samples used for the microarray studies. Figure 4B shows that the expected anatomical registration was confirmed, and that, in addition, there was a significant decrease of 66% in TH mRNA in the MA brains compared with the normal brains (Fig. 1D), consistent with previously published results (Sonsalla et al. 1996).

One of the genes within cluster 2 (Figs. 3D and 4A), the neurofilament light chain gene, Nfl (Yaworsky et al. 1997), was present on the microarrays as two independent spots (AI385738, gene 39; AA253725, gene 50), and this gave an opportunity to assess within-array reproducibility (Fig. 4C). There was a highly statistically significant correlation coefficient for the Nfl expression profiles independently obtained from the two spots on the microarrays, for both the normal and MA brains, confirming excellent within-array reproducibility.

Interestingly, one of the genes in cluster 2, the precerebellin-1 gene (gene 5), is known to be strongly expressed in the cerebellum (Kavety and Morgan 1998). A graphic representation of the spatial expression for precerebellin-1, in both the normal and PD brains, is shown in Figure 4D. The images have strong bilateral symmetry, again emphasizing the good replicability of the technology.

There were some interesting biological relationships within the coregulated clusters of genes. As mentioned above, the neurofilament light chain gene, Nfl, was present in cluster 2 of both the normal and MA brains. Two additional neurofilament genes, neurofilament medium polypeptide, Nfm (W64752; Myers et al. 1987), and -internexin, Ina (AA218283; Chan and Chiu 1996), showed a high level of correlation within cluster 2 of the normal brain, but this relationship was not maintained within the MA brain. Consequently these genes are not shown in Figures 3D and 4A. The pattern of coregulation for the three neurofilament genes within the normal brain presumably reflects their related functions. The relaxation of this coregulation within the MA brain is interesting in light of the fact that aggregated neurofilament subunits are a major protein component of Lewy bodies (LBs), intracytoplasmic inclusion bodies that feature prominently in subcortical neurons of patients with Parkinson's disease (Trojanowski et al. 1998). A bioinformatics analysis found regulatory regions shared between the neurofilament genes, which is shown for two of them (Nfl and Ina) in Figure 5. Also found within cluster 2 are the genes for ring-finger protein 13, Rnf13 (AA189868, gene 38) and ubiquitin-specific protease 9, Usp9X (AA178383, gene 37), both of which are involved in ubiquitin-mediated protein degradation, with ring-finger proteins providing specificity to ubiquitin conjugation (Joazeiro and Weissman 2000).

View larger version (30K): [in this window] [in a new window]   Figure 5   Putative regulatory elements shared between groups of correlated and anticorrelated genes. There were three groups of correlated (+) genes: (1) SIAT9, HIVEP2, SEPP1, SCYD1, NEFL, INA; (2) CUGBP2, MBP, CBLN1; (3) CBLN1, USP9X, PNUTL2; and one group of anticorrelated () genes: (1) MBP, NEFL, CBLN1, PNUTL2, USP9X. The groups are indicated by square brackets. Genes are indicated by UniGene symbol or name (http://www.ncbi.nlm.nih.gov/UniGene). Exons are indicated by short vertical lines and the direction of transcription by horizontal arrows. Homology searches were performed on human homologs of the relevant mouse genes. In all cases, the human gene has the same symbol as the mouse, except for ahuman NEFL = mouse Nfl, and bhuman PNUTL2 = mouse Sept4. The regulatory sequences responsible for correlated expression are shown as squares, those responsible for anticorrelated expression are shown as diamonds. Lines delineate the relationships between the conserved regulatory sequences. Multiple control regions frequently connected the genes. Potential binding sites are: (1) LMO2COM, (2) OCT1, (3) GATA2 and GATA3, (4) MYOD, (5) LMO2COM, (6) NFAT, (7) GFI1, (8) TCF11, (9) HFH2, (10) GFI1, (11) NFAT, (12) FREAC2, (13) OCT1, (14,15) GKLF, (16) TATA, (17) SRY, (18) TH1E47, (19) SRY, (20) HFH2, (21) BARBIE, (22) HFH2, (23) HFH3, (24) HNF3B, (25) HFH2, (26,27) HFH8 and HFH3, (28) HFH-3, (29) GATA, (30) CEBPB, (31) MYOD, (32) LMO2COM, (33) E47, (34) SRY, (35) NFAT, (36) CETS1P54. (*) No known transcription factor binding site.

Gene Networks

To find control regions shared between the correlated and anti-correlated genes of the clusters shown in Figures 3D and 4A, a bioinformatics analysis was performed to look for conserved noncoding sequences (Fig. 5; Table 2). Because of the greater amounts of human genome sequence compared with the mouse, human orthologs of the relevant mouse genes were analyzed. Consequently, only known genes were investigated, and novel genes were omitted from the analysis. Gene pairs were analyzed with gene expression correlation coefficients > 0.75 or <0.75, and the homology search was restricted to sequences 20 kb upstream, 20 kb downstream, and in all introns of the pertinent genes. BLAST was used to find homologies, but not provide reliable appraisal of their statistical significance, because the algorithm uses asymptotic statistical approximations, which are not accurate for shorter sequences (Benson et al. 2000). The resulting homology regions were additionally examined for transcription-factor-binding sites using the TRANSFAC database (Wingender et al. 2000).

The analysis unveiled an intricate arrangement of potential control elements shared between genes, which may play a role in their expression pattern relationships. In all cases, except for three distinct pairs (marked by asterisks in Fig. 5 and Table 2), the potential control sequences harbored known transcription-factor-binding sites. We anticipate that the exceptions may in the future be identified as binding sites for yet uncharacterized transcription factors. Nevertheless, the significance of the potential control regions must still be confirmed experimentally. One striking example that emerged from the analysis is the Nfl gene (human homolog NEFL), which had 14 putative control regions, representing binding sites for no fewer than 19 transcription factors. Interestingly, promoter sequences in the 1.7-kb region upstream of the Nfl start codon have been analyzed using reporter constructs in transgenic mice (Yaworsky et al. 1997). These studies found that the 1.7-kb region contained sequences that allowed initiation of neuronal and myogenic Nfl gene expression in embryos, but not maintenance after birth. Because the analysis reported here examines expression in adult tissues, the uncovered putative control regions may be involved in the maintenance of Nfl expression in the adult brain.

Expression Differences between Normal and MA Brains

In addition to analyses of spatial gene expression in the normal and MA brains, significant (p < 10⁹) global expression differences when averaged across the 40 voxels of each brain were sought for known genes. Scatterplots were constructed that compared the expression levels for these genes between the left and right halves of the normal brain and the left and right halves of the PD brain (Fig. 6A). This analysis showed excellent replicability between the left and right halves for both the normal and diseased brains. Scatterplots were also constructed that compared the left half of the normal brain with the left half of the PD brain, and the right half of the normal brain with the right half of the PD brain (Fig. 6A). These scatterplots revealed that there was about a fourfold difference in expression between the genes repressed or induced in the MA brain. Figure 6B shows the expression level differences for the significantly differentially expressed genes. The results for the left and right halves of the brain are again shown separately and show excellent replicability. There was an approximately equal number of genes induced in the MA brain (16/36) compared with those repressed (20/36), and this did not represent a significant deviation from random (² = 0.22, df = 1, p = 0.64). The differential expression of a number of genes was confirmed using real-time QRT-PCR analysis of the same RNA samples used for the microarray studies (Fig. 6C).

View larger version (37K): [in this window] [in a new window]   Figure 6   Differentially expressed genes. (A) Scatterplots show mean expression levels across the 40 voxels in the normal and MA brains on a logarithmic scale (log2) for genes where p < 109. The comparisons within brains (left and right normal; left and right MA) provide another assessment of replicability (for normal brain, r = 0.90, F[1,34] = 143.98, p < 0.0001; for MA brain, r = 0.98, F[1,34] = 817.76, p < 0.0001). The comparison between brains (left normal and left MA, right normal and right MA) shows the differences in expression between genes repressed and induced in the MA brain. Red line: twofold above equivalent expression level; green line: equivalent expression level; blue line: twofold below equivalent expression level. (B) Mean expression levels across the 40 voxels of the normal and MA brains on a logarithmic scale (log2 ±SEM). (Red) Normal; (blue) MA. The genes are ranked from most (gene 1) to least significant (gene 36; p < 109). To allow replicability to be assessed, the data for the left and right halves of each brain were separated out as the left and right data points, respectively, for each gene. No significant differences were found comparing left normal with right normal (two-tailed t test on log2-transformed data, t = 0.012, df = 70, p = 0.99) and left MA with right MA (two-tailed t test on log2-transformed data, t = 0.33, df = 70, p = 0.75). (C) Confirmation of gene expression differences for selected genes by real-time QRT-PCR analysis of voxel RNA. Expression levels are shown relative to controls and normalized using 18S RNA. The bars indicate mean ± SEM. (***) p < 0.001, (**) p < 0.01, two-tailed t-test. (D) Gene categories in the entire 9000-gene data set (9K) and in those genes significantly different (diff) between the control and MA brains (p < 0.001) when averaged across all 40 voxels. Genes involved in cell/cell interactions (cytoskeleton, extracellular matrix, and cell adhesion) were significantly overrepresented in the regulated subset (p = 0.02), but apoptosis-related genes were not (p = 0.39). The bars show mean ± SD, (*) p < 0.05, as judged using Monte Carlo statistics. The control category represents 15 randomly chosen genes to show the validity of the Monte Carlo analysis. As expected, this category showed no significant difference in frequency between the 9K and regulated data sets (p = 0.48).

Several intriguing genes were revealed by the analysis of expression differences between the normal and MA brains (Table 3), involved in such diverse areas as transcription (e.g., Hdac5, Mxi1, Sp1), cell morphology (e.g., Cdc42, Pkcq), extracellular matrix (e.g., Eln, Lamc2), and signal transduction (e.g., Grb2, Ppp2ca). Interestingly, there appeared to be a relatively large number of genes involved in modulating cell morphology (five genes), suggesting substantial regulation of neuronal outgrowth in response to the loss of dopaminergic neurons in the PD model. Also, many genes were found to play a role in apoptosis (five genes, including Stk2, Ppp2ca, Parg, Siah1a, and qk). This category is not highly ranked in Table 3, as many of the apoptosis-related genes have more than one function and were only assigned to one category in the table.

To more rigorously investigate the classes of genes regulated as a result of the MA treatment, the frequency of genes in the categories described above (cell morphology, extracellular matrix, cell adhesion, and apoptosis) were compared between the entire 9000-gene data set and those genes judged significantly differentially expressed between the control and MA brains (p < 0.001) when averaged across all 40 voxels (Fig. 6D). Using Monte Carlo simulations, genes involved in cell/cell interactions (cell morphology, extracellular matrix, and cell adhesion) were found to be significantly overrepresented in the regulated subset compared with the entire 9000-gene subset. Perhaps surprisingly, apoptosis-related genes were found at an equal frequency in both sets of genes, possibly reflecting the uncertain role of apoptosis in PD (Kösel et al. 1997; Banati et al. 1998; Wüllner et al. 1999). Overall, these data emphasize the importance of cell/cell interactions in response to the neuronal loss of PD.

Other gene categories were also found to be significantly regulated (p < 10^-9) in the PD brains compared with controls (Fig. 6B; Table 3). Sulfation is a major pathway in the biotransformation of many drugs (Xu et al. 2000), and this may explain the induction of the Papss2 (3'-phosphoadenosine 5'-phosphosulfate synthase 2) gene in the MA brain. Interestingly, the Pon2 gene (paraoxonase2, W98586) was highly significantly induced in the MA brain (p < 10⁸; data not shown), and this gene family is also involved in detoxification, as well as being linked with PD (Akhmedova et al. 2001).

The Mtapt/Mapt gene (microtubule-associated protein tau) is repressed in the MA brain. This observation may be related to the fact that tau is mutated in a genetic disorder that features parkinsonism (frontotemporal dementia with parkinsonism on chromosome 17 or FTDP-17), as well as being the constituent of neurofibrillary tangles, the defining cytological lesion for many neurodegenerative movement disorders (Spillantini et al. 2000). Perhaps related to the decreased tau transcript levels, the Ppp2ca gene (protein phosphatase 2a, catalytic subunit,  isoform) is also repressed in the MA brain. Protein phosphatase 2A (Pp2a) binds tightly to tau, and all tau mutations resulting in FTDP-17 decrease the affinity of this association (Goedert et al. 2000).

The qk (quaking) gene (Cox et al. 1999) is repressed in the MA brain and encodes a potential RNA binding protein. Interestingly, mutations in qk result in altered brain dopaminergic signaling (Nikulina et al. 1995), the neurotransmitter system most profoundly affected in PD. Another gene with a connection to dopaminergic signaling is the Sp1 transcription factor, which is repressed in the MA brain and regulates transcription of dopamine receptor genes (Yajima et al. 1998; Hwang et al. 2001).

Singular Value Decomposition Reveals Global Shifts of Gene Expression between Normal and MA Brains

Singular value decomposition (SVD) is a powerful method for parsimonious explanations of complex data sets (Hendler and Shrager 1994; Frackowiak et al. 1997; Alter et al. 2000). This statistical method reduces dimensionality while keeping the maximum possible fraction of the variance from the original data. For example, when used in biomedical imaging, SVD analysis commonly explains data sets in terms of known functional and anatomical boundaries (e.g., cortical vs. subcortical). In the area of gene expression patterns, it might be expected that SVD would show which orthogonal sets of genes (vectors) account for the major variations between the voxels. Essentially, the gene vectors would represent votes for the molecular attributes of the various brain regions in which they are manifest. SVD does not rely on preconceived ideas or hypotheses, and is entirely data-driven.

To ascertain whether SVD would illuminate the large amounts of data from the voxelation studies of the normal and MA brains, we performed an analysis on the anterior/posterior differentially expressed subset of genes. The results are shown in Figure 7, where the first four principal components (PCs), or gene vectors, resulting from the SVD are shown as pseudocolor images. These vector images, rather than representing the expression pattern of any one gene, instead represent multiple (hundreds of) genes, which may be particularly important in regional specification of the brain. The striking bilateral symmetry of these images provides vivid documentation of the excellent reproducibility of the voxelation procedure from an overall survey viewpoint.

View larger version (63K): [in this window] [in a new window]   Figure 7   SVD delineates anatomical regions of the brain. The SVD was performed using the anterior/posterior differentially expressed subset of genes. The spatial patterns resulting from the first, second, third, and fourth PCs are shown. Alongside are the first 400 members of the corresponding gene vectors. The ordinate represents the contribution by the relevant gene to the variation of the vector spatial pattern, and the abscissa represents the genes in decreasing order of differential expression. The pseudocolor scales give the level of expression of the relevant gene vector in the spatial patterns.

The first PC is conserved in the normal and MA brains, and is most strongly expressed in the olfactory lobes and striatum. Similarly, the second PC is conserved between the two brains, and is most strongly expressed in the cerebellum. The third PC, however, is not conserved between the two brains. In the normal brain, this PC is most strongly expressed in the striatum and cerebellum, suggesting a hitherto unsuspected genetic connection between these otherwise functionally and anatomically distinct regions of the brain. In contrast, the third PC of the MA brain shows a dramatic shift away from the striatum and cerebellum toward the hippocampus. This striking alteration in the expression of the corresponding gene vector may be related to the fact that the striatum is the region of the brain most prominently affected by PD. The fourth PC is again conserved between the normal and MA brains, and is strongly expressed in the hippocampus. Together, these findings suggest that the known physiological changes in the striatum of PD brains is accompanied by prominent changes in the corresponding genetic networks responsible for maintenance of this structure.

The reliability of the SVD was assessed by repeating the analysis using a different subset of the data, consisting of those genes with a correlation coefficient between the left and right halves of the brain >0.7, and whose variance of expression through the entire brain was greater than the median of the variance for all genes. The resulting images were essentially identical to those shown in Figure 7, confirming the robustness of the SVD analysis. The bilateral symmetry of the PCs, their robustness, and their restriction to relevant anatomical regions is remarkable, especially considering the twofold uncertainty in the microarray data and the relatively crude spatial maps (40 voxels). With more comprehensive gene surveys and increased resolution, voxelation may ultimately reveal the molecular genetic relationships between various regions of the brain, as well as identifying which areas are most affected in disease. The assignation of gene expression changes to defined neuroanatomical loci may be particularly valuable in identifying responsible brain regions for neuropsychiatric disorders such as schizophrenia (Mirnics et al. 2000), where the location of such regions is presently uncertain.

	DISCUSSION

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS WEB SITE REFERENCES REFERENCES

The investigations reported here show that by combining high-throughput gene expression methodologies with spatial information, voxelation can provide valuable insights not easily obtainable from studies of tissue culture systems. In one example, voxelation helped define networks of anatomically coregulated genes and their relevant control regions. The overwhelming majority of the putative control regions contained potential transcription-factor-binding sites, lending credibility to the analysis. However, transcriptional control regions have considerable degeneracy, and study of genomic sequence from additional organisms may further bolster the significance of the uncovered sequences. The voxelation analysis also allowed identification of genes globally differentially expressed between normal and PD brains, whereas SVD gave insights into the principal neuroanatomical changes of PD at the molecular genetic rather than the traditional histological level.

Significant bilateral symmetry was found in the data using a number of different analytical approaches, suggesting good reproducibility of the voxelation technology. However, the symmetry may also partly reflect the relatively crude spatial maps used for the voxelation, and asymmetry might become apparent at higher resolution in certain situations, for example, as a result of sexual dimorphism (Dluzen and Kreutzberg 1996; Tabibnia et al. 1999). Several studies have used microarrays and targeted dissections to investigate gene expression in the mouse brain (Lockhart and Barlow 2001). Interestingly, in these investigations the cerebellum appeared to be the most distinct brain region in terms of gene expression (Sandberg et al. 2000), reminiscent of the finding in the present study of a posteriorly expressed cluster of coregulated genes.

It should be emphasized that some of the analytical approaches used in this study do not necessarily depend on precisely matched samples. For example, conserved networks of coregulated genes maintained across voxels and between specimens can be identified using human samples (Brown et al. 2002a), where there is an inevitable lack of exact controls. However, analysis of human brains has the decided advantage of disease validity, which may be particularly useful for entities like schizophrenia or bipolar disorder, where the etiology of the disease is obscure and the relevance of the available mouse models is in doubt. In contrast, voxelation of the mouse has the advantages of carefully controlled experiments and matched samples. In the long run, perhaps the most robust understanding will be provided by the judicious combined use of human specimens and model organisms.

Information recovery from voxelation is in principle limited by voxel inhomogeneity and the performance of the analytic tools used to investigate the voxels. Voxel inhomogeneity is a consequence of noninfinitesimal voxel size and hence finite spatial resolution, resulting in decreased signal-to-noise ratios for gene expression. However, this is a ubiquitous phenomenon in brain-imaging technologies (CT, PET, fMRI), but does not prevent these methods from providing useful information. Similarly, useful insights were obtained from the present study, despite the relatively crude spatial maps. The signal-to-noise ratio obtained from finite voxels might be improved by preselecting cells using a lineage-specific marker, such as green fluorescent protein (Peterson 2002). Analytic performance of the tools used to analyze the voxels is limited by sensitivity (qualitative performance: the ability to detect presence or absence of expression), accuracy (quantitative performance: the ability to reliably detect differences in gene expression), and throughput. Microarrays have moderate sensitivity and can only reliably discriminate between twofold differences in gene expression, but have excellent throughput.

Despite information losses in voxelation caused by voxel heterogeneity and analytic limitations, the present study shows that these deficits are overshadowed by the high throughput of the method, which allows for much greater net information recovery than is practicable with classical approaches such as in situ hybridization. This makes the use of voxelation feasible for large-scale study of gene expression in multiple disorders and models, a daunting prospect for in situ hybridization. A further advantage of voxelation is its modality independence, which will allow its use for 3D mapping of the proteome and perhaps even electrophysiology, in addition to investigation of the transcriptome. Nevertheless, an important future goal for voxelation will be to expand the amount of information it provides, both by improvements in resolution and analytic performance.

Although microarrays are a relatively cheap tool on a per gene basis, voxelation will become increasingly expensive as greater numbers of voxels are analyzed in the search for improved resolution in a wide variety of experimental circumstances. It will therefore be important to find ways to drive down costs. Using presently available technology, microarrays are sufficiently sensitive to allow construction of a 325,000-voxel human brain map (Brown et al. 2002a; Peterson 2002), but for such maps to be feasible, the cost per voxel would have to drop by ~25-fold (Peterson 2002). For the mouse brain, microarray technology is sensitive enough to allow construction of ~120 voxel maps, and higher resolution will require improvements in sensitivity. Real-time PCR is sufficiently sensitive to permit construction of 6000-voxel maps of the mouse brain, while still affording opportunities for scalability and automation.

All the goals of improved spatial resolution, analytic performance, and cost will provide substantial challenges. Ultimately, however, high-resolution study of human neuropsychiatric disorders, mouse models of these disorders, as well as the ability to voxelate whole model organisms, will give better comprehension of the logic of the genome, how this logic goes awry in disease, and important starting points for novel therapies.

	METHODS

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS WEB SITE REFERENCES REFERENCES

MA Model of PD

Adult C57BL/6J male mice (10-24 wk, 25-31 g) received four i.p. injections of MA hydrochloride (10 mg/kg per injection using 1.5 mg/mL solution) at 2-h intervals (Sonsalla et al. 1996; Melega et al. 1997). The mice were analyzed 7 d after MA treatment.

Monoamines and Metabolites in Mouse Striatum

Brain samples were weighed wet, sonicated with 250 µL of 0.1 M perchloric acid and centrifuged at 14,000 rpm for 15 min (4°C). The supernatant was filtered through a 0.2-µm PTFE filter and an aliquot diluted with water (1:2) for HPLC analysis. The solid pellet was suspended in 1.0 mL of 0.2 M NaOH for protein assay. For HPLC we used an ESA HPLC Model 580 solvent delivery module (dual-piston pump), and an ESA Coulochem II electrochemical detector with an analytical cell operating at +350 mV and 500 nA. The mobile phase consisted of acetonitrile:sodium phosphate monobasic buffer (75 mM sodium phosphate, 1.8 mM 1-octanesulfonic acid sodium salt [OSA], 12 µM ethylenediaminetetraacetic acid, disodium salt dihydrate [EDTA]) at 9.5:90.5 v/v, pH 3.1 (aqueous phase). The guard column was Adsorbosphere HS, C18, 7.5 × 4.6 mm, 5 µm, and the analytical column was Adsorbosphere HS, C18, 100 × 4.6 mm, 3 µm. Both columns were from Alltech Associates, Inc. The flow rate was 0.8 mL/min, and the injector loop volume was 20 µL.

Voxelation Procedure

Mouse brains were voxelated by first cutting along a transverse plane that included the interaural line and the anterior olfactory lobe, thus producing dorsal and ventral halves. The two halves were placed in a commercially available cutting device (ASI Instruments, RBM-2000C) that incorporates a mold to steady the brain. Machined into the mold are 14 slots separated by 1 mm for precise cutting. In the following description, slot 1 of the mold is most anterior, and slot 14 is most posterior. The brain was placed into the mold such that slice 4 of the mold (corresponding to slice C, Fig. 2) was aligned with the anterior part of the Circle of Willis. Slice 8 of the mold (corresponding to slice G, Fig. 2) was aligned with the point where the anteriormost part of the brain stem/pons could be observed emerging from the ventral surface of the brain. The first and last four slots of the mold were omitted from the voxelation scheme, so that the brain was cut into 10 approximately equal coronal sections of 1 mm thickness, while keeping the dorsal/ventral halves in register. When the blade of the device was removed after cutting each coronal section, the sections, divided into dorsal and ventral halves, remained adherent to the blade, in register, owing to surface tension. At that point each coronal section was bisected down the midline, such that each section consisted of four quarters (superior left, superior right, inferior left, inferior right).

The numbering scheme for the voxels is shown in Figure 2. Slice A was designated as most anterior and slice J the most posterior. Superior right voxels were numbered 1, superior left voxels were 2, inferior right voxels were 3, and inferior left voxels were 4. For example, the inferior right voxel from slice C would be designated voxel C3. Registration with the Mouse Brain Library (Rosen et al. 2000; Williams 2000; http://www.nervenet.org/MBL/mbl.html) was achieved by comparison of sections. Using the inbred C57BL/6J mouse strain, the cutting was found to be essentially invariant from one animal to the next. As judged using the coordinates provided by the Mouse Brain Library, the distance between sections was (mean ± SEM) 1.04 mm ± 0.17 (bregma origin) and 1.04 mm ± 0.18 (interaural origin). To provide sufficient RNA (100 µg) for microarray analysis, equivalent voxels of multiple brains were pooled. Even for the voxels with the smallest amount of tissue, a total of 29 brains provided sufficient material. The smallest voxels were the most anterior (A1-A4), encompassing the olfactory bulbs.

Microarray Analysis

For each voxel, 100 µg of Cy3-labeled voxel RNA and 100 µg of Cy5-labeled control RNA were cohybridized to a separate 9000-gene microarray, as described (Eisen and Brown 1999). The control RNA was used to facilitate interarray comparisons, and consisted of total normal C57BL/6J mouse brain RNA. For each gene, the signal-to-noise ratio was 2.5-fold above background for both the Cy3 and Cy5 channels, and the mean (±SD) signal-to-noise ratio was 7.40 ± 3.25 for Cy3 and 11.04 ± 3.96 for Cy5. The microarray data were processed using two types of normalization procedures. First, spatial trends existing in the data from chip printing were removed by nonlinear transformation of the data sets. The second normalization procedure was designed to compensate for differences in the labeling and chemical properties of the Cy3 and Cy5 dyes, by aligning the histograms of the dye signals both within, as well as between, chips.

Correlation Matrix Clustering

For Figure 3, B and D, the genes in the control omnibus correlation matrix were ordered using a similarity metric. The first row of the matrix was chosen to show a strong contrast between the highest and lowest correlation coefficient for that row. This row was denoted as the base vector, B, with respect to which the remaining rows, R, were arranged in order of decreasing similarity, using a metric consisting of _i(B_i  R_i)², where i is the elements of the rows. Once the omnibus matrix for the normal brain was created, the omnibus matrix for the MA brain was created following the same order, and the data for each brain then separated into left and right halves.

Monte Carlo simulations

The Monte Carlo simulations to assess the similarity of correlation matrices for the normal left with the normal right and the MA left with the MA right (Fig. 3B) used random permutation of the columns of the matrices. The discrepancy between randomly selected pairs of permuted matrices was quantitated using the Frobenius norm of the matrix obtained by subtracting one matrix from the other. The difference between the mean of the resulting distribution and the Frobenius norm obtained from the left and right matrices was then used to assess significance.

Genes in each of the categories shown in Figure 6D were chosen using appropriate keywords (e.g., cytoskeleton, extracellular matrix, and cell adhesion for cell/cell interactions; cell death and apoptosis for apoptosis.) To assess the significance of the gene frequencies, Monte Carlo simulations were performed using two data sets: the entire 9000-gene data set and the subset consisting of those genes differentially expressed (p < 0.001) when averaged across all 40 voxels of the control and MA brains. The simulations were performed by randomly discarding the data for 30% of the genes in each data set and recalculating the mean. Significance values were assigned for each gene category by calculating the area of overlap for the two distributions from each data set. To show the validity of the Monte Carlo simulations, a control category was analyzed, consisting of 15 genes chosen randomly that were common to both data sets. As expected, no significant difference was found between the data sets for this category (Fig. 6D).

Real-Time QRT-PCR

Real-time QRT-PCR was performed using TaqMan (Gibson et al. 1996) and One-Step RT-PCR Master Mix following the manufacturer's instructions (Applied Biosystems). Reverse transcription used 100 ng of total voxel RNA for TH, and 10 ng for Ecm1, Ap1b1, Psme1, S100a6, and Stk2. Detection used an Applied Biosystems Prism 7700 Sequence Detector (Perkin-Elmer), and the data were analyzed using SEQUENCE DETECTOR software. Normalization used 18S RNA. Genomic contamination was excluded (Brown et al. 2002b) by the use of primers that cross the intron of the housekeeping gene GdX (Filippi et al. 1990) and no reverse transcriptase controls.

Images

For the anatomical images of gene expression (Fig. 4D) and SVD (Fig. 7), a midlevel transverse section was used (5.40 mm interaural, 4.60 mm bregma, section 8 of the Mouse Brain Library), and expression levels were averaged between corresponding dorsal and ventral voxels on the same side of the brain. All the images used a data set consisting of those genes that were most strongly differentially expressed between the anterior half (20 voxels) and posterior half (20 voxels) of the brain (p < 0.05). Because the cerebral cortex is featured in approximately equal amounts in both the anterior and posterior halves, it was expected a priori that the images would principally show features in parasagittal locations rather than in the periphery, where the cortex is located. Consequently, a Bayesian approach to image creation was used, based on a kernel that gave preference to the center of each voxel, and thus used prior probabilities from the known anatomy of the voxelated brain. Furthermore, a prior assumption of nearest-neighbor continuity was used, resulting in smoothed expression patterns over the voxels.

	WEB SITE REFERENCES

TOP ABSTRACT INTRODUCTION RESULTS DISCUSSION METHODS WEB SITE REFERENCES REFERENCES

http://www.nervenet.org/MBL/mbl.html; Mouse Brain Library.

http://www.ncbi.nlm.nih.gov/geo; Microarray data under the series accession no. GSE30.

http://www.ncbi.nlm.nih.gov/UniGene; UniGene.

http://www.pharmacology.ucla.edu/smithlab/genome_multiplex; All study results and supplementary material.