Many methods have been developed to standardize the rodent brain map (e.g., [27], [32], [33], [35], [41], [42]). In particular, Gabbott et al. (2005) reported a method similar to the cortical box method to investigate the cortical projections from the rat frontal areas [35]. Our method is also conceptually similar to the surface- based atlases developed by Van Essen and coworkers [43]. The advantage of these methods is that by flattening the three- dimensional cortex, it becomes easy to understand intuitively the global picture of gene expression. It is also important that the cortical box method enables the quantitative analyses of the spatial distribution data. Even when there are sample-to-sample variances owing to various reasons, we can extract useful information by averaging the data. We can also estimate the significance of such variances (see Fig. 6B). We noted that the variability of gene expression were concentrated at area borders (Fig. 6B), which could be attributable to individual difference of the area architecture, although it is possible that such variance was derived from experimental variance. One advantage of our method over other flattening techniques is the simplicity and ease of use, because it is optimized for the simple sulcus structure of the rodent posterior cortex. On the other hand, it would be difficult to directly apply it to the convoluted cortex of other mammalian species. For example, because the thickness of each lamina varies greatly in areas for primate brains (e.g., see Fig. S3 of [23]), it will become difficult to construct layer maps using the same layer borders for different areas. However, if we are to limit the analysis to a subregion of the cortex that can be defined by clear-cut landmarks (such as sulci), cortical box method may be used to analyze convoluted brains as well. Despite the many useful features of the cortical box method as pointed above, caution is required when interpreting the data, because it does not offer information on the expression at the cellular resolution level. For example, in the layer 6 maps of ER81 and Nurr1, the maps of two genes overlapped in the laterocaudal areas (Figs. 5 and 6, TeV and Ect), although the double ISH results demonstrated that they were not coexpressed (Fig. 3G–I). It is also notable that, owing to normalization, the low expression levels tend to be ignored in a global view. For example, although RORbeta mRNA is clearly expressed by layer 5 neurons throughout areas (Fig. 2), it does not show up in the layer 5 map (Fig. 6A). In the case of ER81, the very high expression level in the medial areas obscures its widespread distribution across the entire neocortical areas in layer 5. The cortical box method needs to be coupled with careful analyses of expressions at the cellular level. It should also be noted that, by fitting the cortex into rectangle, the information of the cortical thickness is lost, although the relative expression value is preserved throughout layers. To examine the differential cortical thickness across areas, other analytical method needs to be used. Identification of cortical areas by cortical box method In previous studies, cortical areas were successfully delineated in an objective manner by Nissl-GLI analyses [8], [36–38]. Therefore, we relied on Nissl-GLI patterns to identify cytoarchi- tectonic areas in our standard three-dimensional space. The layer 4 map of Nissl-GLI was consistent overall with the Nissl-GLI map of Zilles and coworkers [8]. However, there were several points wherein we incorporated the area classification by other researchers. For example, we found a subarea ‘‘Par1L’’ in the lateral region of the primary somatosensory cortex, which is not shown in the original map of Palomero-Gallagher and Zilles [8]. Judged from the location, Par1L seems to correspond to the representation region of the upper lip, which was recently noted by others [39]. We also found another subarea ‘‘dysgranular zone ISH Profiling of Rat Cortex (DZ)’’ between areas HL and Par1. This subarea appears to correspond to the most medial zone of a matrix of dysgranular cortex into which the barrels are embedded [44–47]. The DZ appears to partially overlap with ‘‘FL’’ in the map of Palomero- Gallagher and Zilles [8]. Although the borders of DZ were not determined by GLI analysis [8], our data showed that they extended anteroposteriorly along HL. The consistency with the RORbeta gene expression supports this area delineation. In the auditory cortex, we distinguished Te1 and Te3V, on the basis of Nissl-GLI and RORbeta expression. Te3V in our map is considered to correspond to the belt region of the auditory cortex [8], [48]. Ventral to TeV, we delineated Ect according to Paxinos and Watson [49] and Swanson [50]. We emphasize that the cytoarchitectonic areas determined by Nissl-GLI are well consis- tent with the expression patterns of the layer-specific genes, validating the area demarcation by our method. Significance of area-specific gene expression Previous studies using receptor autoradiography have revealed that the brain’s chemoarchitectonic organization is correlated with cyto- and myeloarchitectonical organizations [7], [51]. Our analysis also revealed tight correlation between gene expression and the cytoarchitectonic area. In particular, RORbeta and ER81 patterns were quite similar to those of Nissl-GLI in layers 4 and 5, respectively. This is anticipated to some extent, because both Nissl- GLI and gene expression intensity should positively correlate with the neuronal density. However, RORbeta and ER81 mRNAs were expressed by a subpopulation of cortical neurons and the variations in their densities appeared to be much greater than those expected from the neuronal density of the layer of interest (Fig. 2). Besides, the perceived intensity of labeling per cell also varied across areas. Regarding Nurr1 gene expression, the positive cells represent only a minor population of the layer 6 neurons, and the area distribution pattern was completely different from that of the Nissl-GLI. These observations suggest that there may exist some rules that commonly affect their area-specific expression patterns other than neuronal density. By using PCA, we tried to find such rules and obtained two principal components (PC1 and PC2). The PC1–PC2 plot of different cortical areas shown in Fig. 7B indicates the features of these areas characterized by RORbeta, ER81 and Nurr1 gene expressions. This figure clearly shows that the primary sensory areas, Par1, Par1L, Te1, Oc1M and Oc1B, cluster together at lower PC1 values, while the higher sensory areas and multimodal areas are away from the ‘‘cluster’’ of primary sensory areas with higher PC1 values, but with various PC2 values. This pattern demonstrates the similarity of gene expression in the primary sensory areas, as well as the diversity of gene expression patterns in the higher sensory areas. There are several possible explanations for the similarity of the primary sensory areas. First, the primary sensory areas generally have ‘‘granular’’ layer 4, which contains a high density of thalamor- ecipient neurons. The PC1 score may be positively or negatively correlated with the neuronal density as we discussed earlier. Second, the primary sensory areas receive strong inputs from the primary sensory thalamus [52], [53], which can be visualized by cytochrome oxidase staining [54], [55]. It is well known that certain genes exhibit activity-dependent regulation during devel- opment or in the adult [17], [56–58]. The expressions of RORbeta, ER81 and/or Nurr1 genes may be affected, directly or indirectly, by the thalamocortical inputs and contribute to the low PC1 scores in the primary sensory areas. In this context, it is interesting that area DZ exhibits a high PC1 score despite its location within Par1. Previous studies suggest that the dysgranular PLoS ONE | www.plosone.org 9 September 2008 | Volume 3 | Issue 9 | e3266