A diffusion MRI-based spatiotemporal continuum of the embryonic mouse brain for probing gene-neuroanatomy connections.

The embryonic mouse brain undergoes drastic changes in establishing basic anatomical compartments and laying out major axonal connections of the developing brain. Correlating anatomical changes with gene-expression patterns is an essential step toward understanding the mechanisms regulating brain development. Traditionally, this is done in a cross-sectional manner, but the dynamic nature of development calls for probing gene-neuroanatomy interactions in a combined spatiotemporal domain. Here, we present a four-dimensional (4D) spatiotemporal continuum of the embryonic mouse brain from E10.5 to E15.5 reconstructed from diffusion magnetic resonance microscopy (dMRM) data. This study achieved unprecedented high-definition dMRM at 30- to 35-µm isotropic resolution, and together with computational neuroanatomy techniques, we revealed both morphological and microscopic changes in the developing brain. We transformed selected gene-expression data to this continuum and correlated them with the dMRM-based neuroanatomical changes in embryonic brains. Within the continuum, we identified distinct developmental modes comprising regional clusters that shared developmental trajectories and similar gene-expression profiles. Our results demonstrate how this 4D continuum can be used to examine spatiotemporal gene-neuroanatomical interactions by connecting upstream genetic events with anatomical changes that emerge later in development. This approach would be useful for large-scale analysis of the cooperative roles of key genes in shaping the developing brain.

The mammalian brain undergoes rapid development during embryogenesis, which involves neuronal proliferation, migration, differentiation, and the formation of neural circuits and functional units, and these processes are driven by multiple intricate yet highly dynamic processes, each controlled by a large network of genes. Indeed, identifying gene-expression patterns in the brain is the first step to understanding morphogenesis and patterning of the developing brain. Developmental neurobiologists have spent decades mapping the gene-expression patterns in the mouse brain and other model systems, mainly based on cross-sectional in situ hybridization (ISH) data. Recent advances in gene technology and high-throughput ISH, such as those used to create the Allen Developing Mouse Brain Atlas (ADMBA) (1), have facilitated the collection and mapping of three-dimensional (3D) gene-expression data in a common anatomical framework, which permits systematic examinations of the cooperative roles of multiple genes in shaping neuroanatomical phenotypes. The inextricable nature of spatial and temporal patterns in both gene expression and neuroanatomical phenotypes, however, makes it challenging to characterize their dynamic relationship based solely on discrete datasets, as early changes in gene expression may result in anatomical changes only becoming detectable at later stages (2). Instead, both anatomical and gene-expression data should be treated as four-dimensional (4D, spatial and temporal) entities in a combined spatiotemporal domain, or the developmental continuum, so that upstream events that trigger anatomical changes later in development can be linked.

Building the developmental continuum demands new toolsets to capture neuroanatomical changes of the developing brain and to connect them with gene-expression data over time. Advanced neuroimaging tools have been developed to capture developing brain phenotypes, both in vivo and ex vivo (3–6). MRI, with its rich soft tissue contrasts and abilities to characterize intact neuroanatomy in 3D without sectioning and staining, has been used to quantitatively map the structural and functional development of the embryonic mouse brain (6–10) as well as to detect alternations due to genetic and environmental factors (11–13). In particular, diffusion MRI (dMRI) provides superior contrasts in the developing mouse brain compared to more commonly used T1- or T2-weighted anatomical imaging (14). dMRI contrasts not only enable delineation and tracking of rapidly changing structures but also can capture tissue microstructural organization and structural connectivity in the embryonic and fetal brains (15–18). Due to the small size of embryonic mouse brains (∼1.5 mm3 at E10.5), high-resolution dMRI, or diffusion MR microscopy (dMRM), is essential for building the developmental continuum. Advanced computational brain-mapping techniques have been developed to connect MRI data acquired at different developmental stages (19, 20) and to directly relate MRI and gene-expression data, overcoming the vastly different spatial resolution between them ((101 to 102 µm) versus ((10−3 µm)) (21, 22). Several studies have aligned anatomical MRI and ISH images in the Allen Mouse Brain Atlas (AMBA) (23, 24) space to identify the morphological correspondence of genes in the adult mouse brain (25–27). However, the spatiotemporal interaction between gene expression and developing neuroanatomy in the embryonic or neonatal brain remains a largely unexplored territory.

In this study, we reconstructed a 4D continuum of the developing embryonic brain of the CD1 mouse based on dMRM data at a ultra high resolution of 30 to 45 µm isotropic. We quantitatively examined microstructural and morphological changes in the developing continuum, based on which, regional clusters with distinct developmental trajectories were identified using computational anatomy approaches (28, 29). We further demonstrated the feasibility of integrating 4D gene-expression data into the developmental continuum to allow direct visualization of their interactions in both normal and abnormal developing mouse brains.

Results

dMRM Enables Tracking of Gray and White Matter Development in the Embryonic Mouse Brain.

Ex vivo dMRM data were acquired from embryonic mouse brains (CD1, E10.5 to E15.5, and n = 5 at each stage) at isotropic resolutions of 30 μm, 35 μm, and 45 μm for the E10.5, E11.5 to E14.5, and E15.5 embryonic mouse brains, respectively (Fig. 1). The high resolution together with the unique contrasts from dMRM allowed us to characterize and track several key aspects of embryonic brain development. Changes in gross brain morphology can be appreciated from T2-weighted anatomical images (Fig. 1, Upper), which show how the five primitive vesicles at E10.5 expand rapidly during E10.5 to E11.5 to form the basic layout of brain compartments at E12.5. For tracking the development of internal brain structures, dMRM provides sharp gray and white matter contrasts (e.g., fractional anisotropy [FA] maps in Fig. 1) and information on tissue microstructural architecture as shown by tract density images (TDI) in Fig. 1. At the 30-µm resolution, the thin telencephalic wall at E10.5 can be delineated with radial streamlines in the TDI map, reflecting the cellular scaffold formed by radial glial fibers in the neuroepithelium (NE). The development of cortical laminar structure is comprised of the NE and the later-emerging cortical plate (CP) and intermediate zone (IZ), which can be distinguished based on their distinct orientational profiles from the dMRM results (Fig. 1). The sequential development of the IZ that is initiated at the pallio–subpallial boundary of the ventropallial prepiriform area at E12.5 and then extended into the dorsal pallium/isocortex area at E14.5 (), possibly reflects the tangential migration of interneurons from subpallium to the cortical regions (30, 31).

Fig. 1.

dMRM of the developing embryonic mouse brain. (A) Image acquisition. (Upper) The embryonic mouse brain was embedded in agarose and imaged using a saddle coil on an 11.7-Tesla scanner. (Lower) A custom diffusion-weighted gradient and spin-echo sequence was used to achieve 3D dMRM acquisition with a 20-fold acceleration (35). (B) Representative T2-weighted (T2w) and FA maps of embryonic mouse brains acquired using this setup. (C) TDIs visualize microstructural organization in the embryonic mouse brains using streamlines. The color indicates the direction of the streamline: red, medial–lateral; green, dorsal–ventral; and blue, rostral–caudal. (C') Zoom-in illustration of the E13.5 mouse cortex. The NE, IZ, and CP can be clearly separated in the TDI map due to their distinct microstructures as illustrated by the 3D schematics. The NE and CP show a radial organization, whereas the IZ contains both radial and transverse fibers, as indicated by a zoom-in view of the estimated fiber orientation density map. (D) Reconstruction of early white matter tracts. (Upper) Three-dimensional rendering of the embryonic brains from E10.5 to E15.5 and development of thalamus-related fiber bundles from E12.5 to E15.5. (Lower) Development of projection and commissural fibers along the developmental timeline. Abbreviations: stria medullaris: sm, stria terminalis: st, fasciculus retroflexus: fr, cerebral peduncle: cp, optic tract: opt, anterior commissure: ac, fimbria: fi, medial lemniscus: ml, and mammillothalamic tract: mt.

Early fiber tracts reconstructed from dMRM data from as early as E12.5 provide key landmarks for tracking the structural connectivity in the embryonic brain (Fig. 1). The first group of visible tracts included the stria medullaris, stria terminalis, and fasciculus retroflexus (Fig. 1, Upper), which connect the thalamus with the hypothalamus, septum, and amygdala. They are followed by major projection and association tracts, such as the cerebral peduncle at E13.5, the optic tract at E14.5, and the anterior commissure and fimbria at E15.5 (Fig. 1, Lower). Several smaller tracts, such as the medial lemniscus and mammillothalamic tract, emerged later at E15.5. The development of these fiber bundles closely follows the timeline of brain development (32, 33). More-comprehensive views of these fiber tracts are available in Movies S1 through S4. Taken together, dMRM provides rich anatomical information of both gray matter and white matter for characterizing early brain development.

Building a 4D Mouse Embryonic Brain Developmental Continuum.

From individual embryonic mouse brain images at each stage, group-averaged population atlas images were generated using an iterative registration procedure (34, 35) (Fig. 2). The resulting atlas brains () were unbiased representations of average anatomical features at each stage with improved signal-to-noise ratio compared with those of individual brains, while preserving sharp contrasts for structural delineation (Fig. 2). We then computed diffeomorphic mappings between atlas images at adjacent stages from E11.5 to E15.5 (Fig. 2). The diffeomorphic constraint ensured that the large deformations necessary to capture the drastic morphological changes during development would not disrupt the topological relationships among brain structures (36). We excluded the E10.5 brains in the following computational anatomy analysis due to drastic change in brain morphology from E10.5 to E11.5, which makes accurate mapping between these two stages challenging. The remaining subtle topological changes (e.g., newly-emerged axonal tracts) do not disrupt the mappings due to the imposed smoothness condition.

Fig. 2.

Quantification of 4D embryonic brain development using computational anatomy techniques. (A) The computational anatomy pipeline. The dMRI metrics of individual mouse brains (n = 5) were coregistered to generate an average atlas at each embryonic day from E10.5 to E15.5 using iterative registration procedures. Sequential registrations between consecutive embryonic stages were performed to obtain the transformation matrices and GVFs as well as to align all brains in the E13.5 common space for downstream analysis. (B) The diffusion tensor–based, directionally encoded colormaps of individual mouse brains (Upper) and the population-averaged atlases (Lower) from E10.5 to E15.5. (C) The GVF maps of morphological changes between neighboring stages along with zoom-in views of the GVs in the cortical region. The colormaps underlying the GVF indicate the directions of the growth vectors with color codes shown on the Left. Color codes: red, medial–lateral; green, dorsal–ventral; and blue, rostral–caudal.

After spatially normalizing these cross-stage mappings to a common space on an E13.5 brain atlas, the midpoint, we concatenated them into a 4D growth vector field (GVF) (Fig. 2). This massive GVF encodes tissue displacement and deformations over the entire brain from E11.5 to E15.5 and allows us to follow the morphological developmental trajectories of each brain region over time, both forward and backward. For example, in the dorsal region of the cortex (yellow boxes in Fig. 2), the GVF shows the directions of the expansion of the telencephalic wall both within the coronal plane (red–green color codes for E11.5 to E14.5) and along the rostral–audal axis (blue color for E14.5 to E15.5). By interpolating the 4D GVF to a finer time interval (0.1 d), we can visualize the continuous developmental processes ( and Movie S5). Together, the dMRM atlases and the 4D GVF formed the developmental continuum of embryonic mouse brain development.

Serial log-Jacobian determinant (log-Jt) maps calculated from the 4D GVF (Fig. 3) are used to visualize the local tissue deformations (expansions or shrinkages) over time for quantifying brain morphological changes. The cerebral cortex shows the fastest expansion among all brain regions during development, but its spatial pattern of expansion evolves over time. For example, while the entire cortex expands in a relatively uniform fashion from E11.5 to E13.5, active growth primarily localizes to the temporal region during E13.5 through E15.5 (white arrows in Fig. 3). Morphological changes associated with emerging white matter tracts are also present in the map, such as the anterior commissure and the cerebral peduncle (red arrow and black arrows in Fig. 3) showing high log-Jt values from E14.5 to E15.5. Moreover, the primary direction of morphological change and the degree of change along the primary orientation can be quantified by strain tensors (37) computed from the GVF. From E11.5 to E14.5, the primary orientations of strain tensors in the embryonic cortex are mostly parallel to the cortical surface, suggesting that tangential elongation of the cortical surface is dominant here, whereas from E14.5 to E15.5, the cortical tensors become perpendicular to the surface, indicating that thickening of the CP may become dominant during this stage (Fig. 3). The change of tensor orientations can be clearly visualized from the polar plot in Fig. 3 that shows that the tensors are close to 0/180° (tangential to the cortical surface) at E11.5 and gradually move to 90° (perpendicular to the cortical surface) at E15.5.

Fig. 3.

Quantification of morphological brain development based on the extent (log-Jt) and shape (strain tensor) of GVF between consecutive developmental stages. (A) Log-Jt maps of the brain deformation from one stage to the next. The maps are rendered in 2D and 3D views, with the color bar indicating the extent of deformation. Note that the different color bars are used for different stages to emphasize regional contrasts. The red arrow points to the emerging anterior commissure, the black arrows point to the developing cerebral peduncle, and the white arrows point to the temporal cortex that experiences fast growth during E13.5 through E15.5. (B) Strain tensor maps showing the primary orientation of the morphological development from one stage to the next, with zoom-in views of the cortical regions. (C) A polar plot of the primary orientations of strain tensors of morphological development in the embryonic cortex at different stages. Each dot represents a cortical voxel, and its location in the plot shows orientation of the strain tensor in this voxel with respect to the cortical surface (0/180°: tangential and 90°: perpendicular) and magnitude (distance to the center).

Mapping Gene-Expression Data to the Developmental Continuum.

We transformed gene-expression data to the dMRM-based 4D continuum using a pipeline of landmark-based affine and diffeomorphic transformations () followed by intensity-based registrations (Fig. 4). ISH data of 10 different genes were obtained at E11.5, E12.5, E14.5, and E15.5, including Emx2, Wnt7a/7b, Robo1, Slit2, etc., which are known to play key roles in cortical, subcortical, and white matter development in the embryonic mouse brain (38–45). The spatiotemporal patterns of these genes are visualized in the dMRM atlas space (Fig. 4). In order to quantify the spatiotemporal correlation between gene expression and dMRM-based neuroanatomy, we measured regional expression/dMRM profiles by segmenting the embryonic brains into 26 prosomeric regions () according to the ADMBA definition. The segmentation-based analysis showed region-specific changes of microstructural and morphological measurements across E11.5 through E15.5 ().

Fig. 4.

Integration and correlation between gene expression and dMRM data in the developmental continuum. (A) Pipeline for dMRM–ISH registration and subsequent quantitative analysis. Series of ISH slides were transformed to the dMRM-based space at the corresponding stage using landmark-based affine and diffeomorphic registrations, such that multiple genes can be analyzed in a common space. The embryonic brain was then parcellated into prosomeric regions by coregistration with ADMBA, based on which, the spatiotemporal correlation between the gene expression and dMRM data could be analyzed. (B) Three-dimensional views of the gene expression in the dMRM-based space for 10 genes from E11.5 to E15.5. (C) Spatiotemporal patterns of gene expression in 17 forebrain segments (telencephalon and diencephalon) at E11.5, E12.5, E14.5, and E15.5. The sum of expression intensity in each segment was first log transformed and then normalized across the brain into 0 to 1. Note that the midbrain and hindbrain segments were not included as some of the ISH sections did not cover these regions. (D) Spatiotemporal patterns of the dMRM-based FA and MD metrics in the forebrain prosomeric segments across gestation. (E) Spatiotemporal correlations between gene expression and dMRM profiles. The color code indicates the correlation coefficient (r) and * highlights the significant correlations with |r| > 0.2. Abbreviations: Pall: Pallium, Spall: Subpallium, Hy: hypothalamus, and Dienc: Diencephalon.

Two-dimensional spatiotemporal profiles of gene expression in the prosomeric regions across E11.5 through E15.5 showed that some genes (e.g., Emx2 and Slit2) exhibit localized expressions that decrease with age, while others (e.g., Wnt7b and Robo1) have more-widespread distributions that generally increase with age (Fig. 4). Two-dimensional profiles of dMRM metrics showed a decrease of mean diffusivity (MD) across the forebrain from E11.5 to E15.5, likely due to the increasing microstructural complexity and increased restrictions on water diffusion (46, 47), whereas the FA values exhibited both increasing and decreasing trends in different brain regions (Fig. 4). Correlations between these 2D profiles showed a potential link between gene expression and neuroanatomical structures (Fig. 4). We found that spatiotemporal pattern of FA is positively correlated with Emx2 and Wnt7a (r > 0.2 and P < 0.05), possibly relating to the decreasing expression of these genes that are most active at early embryonic stages for neurogenesis and cortical formation (38, 48), while the diffusivity measurements showed noticeable negative correlations with Wnt7b, Robo1, Lhx6, and Pak3 (r < −0.2 and P < 0.05), possibly because the age-dependent reduction in diffusivity is opposite to the increased expression of these genes during axonal outgrowth and neuronal maturation at later stages (39, 43).

Connecting Genetic Events with the Developing Neuroanatomy.

Mapping gene expression in the 4D continuum allows direct examination of expression patterns of single or multiple genes in the context of developing neuroanatomy. For example, the expression patterns of Emx2 closely followed the NE layer of the developing cortex from E11.5 to E15.5 (Fig. 5), consistent with the known role of Emx2 in regulating cortical development and regionalization (38). Wnt7a and Wnt7b, known to regulate synaptic formation and dendritic arborization of hippocampal neurons (39, 40), exhibited complementary spatial localization in the developing brain (Fig. 5), consistent with their roles in regulating dorsal–ventral brain patterning in conjunction with Bmp, Fgf, and Shh signaling (49).

Fig. 5.

Genetic–neuroanatomy interactions in the developmental continuum. (A) Emx2 expression overlaid onto FA maps from E11.5 to E15.5, showing that it is selectively expressed in the pallium of the E11.5 and E12.5 brains (mainly the NE layer) and the lower layer of the E14.5 and E15.5 cortex (NE layer, red arrows). Gene expression intensities were normalized to 0 to 1. (B) Wnt7a (red) and Wnt7b (blue) expression overlaid on FA maps, showing their complementary spatial preferences with Wnt7a primarily in the ventral and medial hippocampus and Wnt7b in dorsal and lateral regions. (C and D) Slit2 (C) and Robo1 (D) expressions transformed from E12.5 (Upper) to E15.5 (Lower). The yellow arrows point to the anterior commissure at the midline. Three-dimensional renderings of the transformed E12.5 gene expression and fiber pathways at E15.5 show their spatial relationship. More-comprehensive views can be found in Movies S6 and S7.)

More importantly, the developmental continuum provides a way to trace earlier gene expressions at the site of anatomical changes that only become visible at later stages by mapping gene-expression profiles of interest along the temporal axis using the GVF. Take the Slit-Robo signaling as an example, which is known for its key role in formation of commissural fibers (50, 51). In the dMRM data, fibers in the anterior commissure were first seen crossing the midline at E15.5. By transporting the Slit2 and Robo1 expressions at E12.5 to the E15.5 mouse brain, we were able to examine earlier expression patterns of Slit2 and Robo1 with the anterior commissure at E15.5 in the same space (Fig. 5 ). The observation that earlier Slit2 expression, flanked by Robo1 expression, aligned with the location where the anterior commissure later crosses the midline is consistent with existing knowledge on the role of Slit2 in channeling axons toward the midline and Robo1 in mediating the postcrossing guidance of commissural fibers (52).

We further demonstrate how the proposed developmental continuum can be used to trace the abnormal brain development in the Sas4 embryonic mouse brain, which shows abnormal progenitor cell distribution (Fig. 6) (53). dMRM revealed altered microstructural organization in the IZ of E15.5 Sas4 embryos compared to wild-type (Fig. 6 ), possibly due to the displaced radial glial progenitors (RGPs) in the mutant brain. By tracing backward in time within the developmental continuum to E12.5, we found that the abnormal region at E15.5 originated from approximately the center of the NE layer at E12.5 (Fig. 6), matching with the Sas4-induced centriole deficits and cell death at this stage (Fig. 6) (53).

Fig. 6.

(A–C) Tractography-based analysis of the Sas-4 mouse brain revealed reduced axons (B) and impaired fiber orientation density (C) running tangential to the cortical surface in the IZ compared to the wild-type mice, possibly due to the displaced RGPs in the mutant brain (A). (D–E) The abnormal loci detected at E15.5 can be traced back to E12.5 based on the GVFs in the developmental continuum, which matched well with the Sas4-induced centriole deficits at E12.5. Note that the p53 mechanism does not take effect until E15.5, and thus, Sas4 knockouts are essential the same as Sas-4 at E12.5. Immunohistology stained sections in A and D are adapted from ref. 53 with permission. (Scale bars, 25 µm in D, Top, and 50 µm in A and other rows of D.)

Defining Regional Clusters with Distinct Developmental Trajectories.

We further used an unbiased data-driven approach to search for regional patterns of brain development by clustering the voxels in the 4D continuum that shared similar temporal trajectories in terms of FA and log-Jt. Six clusters, or intrinsic modes, with distinct developmental trajectories were identified, of which five clusters corresponded to anatomically well-defined regional domains (Fig. 7 ). Voxels in the different clusters are well separated in the principal component feature space (Fig. 7). Mode 1 primarily locates in the NE layer of the cortex and the cerebral peduncle, where the FA shows a biphasic change with a peak at E13.5 and the log-Jt keeps decreasing from E11.5 to E14.5. Mode 2 covers the CP layer of the cortex and several subcortical regions, with similar FA and log-Jt trends as Mode 1 except that the FA values are much lower. Mode 3 includes a large portion of the ventral forebrain, including part of the regions corresponding to the adult hypothalamus and amygdala, with relatively stable, low FA values and rapidly decreasing log-Jt. Mode 4 mainly involves the ganglionic eminence and ventricular zone of the hypothalamus, with FA and log-Jt both decreasing with gestation. Mode 5 is selectively associated with the thalamus, which shows decreasing FA and log-Jt. Mode 6 covers the rest of the brain (e.g., the ventricular areas) and is excluded in the following gene-expression analysis. We varied the number of clusters from four to eight, and a cluster number of six resulted in the most anatomically interpretable outcomes (). The relation between the ADMBA-based prosomeric segmentation and the developmental modes can be found in .

Fig. 7.

Developmental mode analysis to find regions with similar developmental trajectories. (A) K-means clustering of the whole-brain voxels into six modes based on the developmental trajectories of FA and log-J from E11.5 to E14.5. The voxels are plotted in the feature space based on the first and second principal components (PC1 and PC2). (B) Two-dimensional and volumetric rendering of the six anatomically meaningful developmental modes. Note that Mode 6 that covers mostly nonbrain tissues was not included in the following analysis. (C) Individual developmental modes and corresponding developmental trajectories of FA and normalized log-J, plotted as the mean and SD of all voxels within the developmental mode. (D) Two-dimensional dendrograms of the expressions of 10 coregistered genes in each developmental mode from E11.5 to E15.5.

Further examination of gene expression across the five modes suggests that some of the genes selectively resided in certain modes and such preference changed over time (). We quantified the gene-expression profiles in individual modes using 2D dendrograms (Fig. 7), which showed that Modes 1, 2, and 4 shared the same clustering pattern along the gestation age axis. In Mode 1, Wnt7a, Pak3, and Emx2 are clustered together as they are all highly expressed at E11.5 and E12.5, and their expressions drop at later ages, while in Modes 2, 4 and 5, Wnt7a, Pak3, and Lhx5 are closely clustered due to their high expression throughout development.

Integrated Analysis with ADMBA.

Using the anatomic gene expression atlas (AGEA) tool in ADMBA, which groups regions with similar genetic expression across ∼2,000 genes (54) (), we found that three of the five developmental modes (1, 3, and 5) have direct correspondence to regional clusters of gene expression (, Fig. S9 ), indicating the regions with similar anatomical changes are potentially driven by a similar set of genes. We further performed gene enrichment analysis () on the set of genes with similar spatial distributions in Mode 1, 3, and 5, respectively. Results suggest that these genes are mostly related to neurogenesis, differentiation, migration, axonogenesis, and axon guidance (), and moreover, regional specificity could be identified (e.g., the genes enriched in subpallial [Mode 3] are more closely related to axon development compared to those in pallidum [Mode 1]) (, Fig. S9).

Discussion

During the past decade, the development of high-throughput ISH platform and rich genetic resources, such as the ADMBA (1), have advanced our understanding of the genetic mechanisms controlling mammalian brain development. In the meantime, technical advances in MRI have enabled virtual dissection of embryonic neuroanatomy and even imaging of living embryos (10, 55–57). Studies on the relationships between gene expression and imaging-based structural/functional phenotypes are potentially useful for screening candidate genes that play critical roles in normal and genetically modified mouse brains. While several studies in the adult mouse brain have demonstrated remarkable links between neuroanatomy and gene-expression patterns (25–27), the dynamic relationships between gene-expression patterns and anatomical structures in developing mouse brain are not well understood. By adding a temporal dimension to the analysis of structure-gene relationships in the developing brain, we demonstrated that, for the first time, dynamic patterns of gene expression can be integrated and traced forward or backward in time in the context of changing neuroanatomy based on ultra high-resolution dMRM data. This approach not only provides microstructural, tractography, and morphological details of the embryonic brain, but also facilitates the search for genotype-phenotype associations across the 4D continuum, which is important when the anatomical change of interest reflects the cumulative results of multiple genes over time.

dMRM Provides a Mesoscopic View of the Developmental Continuum.

Compared to relaxation-based MRI contrasts, which are sensitive to tissue chemical compositions, contrasts in dMRI hinge on tissue microstructural organization, such as the orientations, densities, and sizes of axons and radial glia. These contrasts are uniquely suited for imaging unmyelinated embryonic brains (9), as several dMRI-based studies have demonstrated in the developing rodent (15, 16, 58), primate (17), and human brains (18, 59). In order to capture emerging brain structures and provide sufficient anatomical details for computational analysis, we used MRI pulse sequences (35) and hardware specially designed for microimaging to achieve an isotropic resolution of 30 µm, which represents a five- to 15-fold reduction in voxel size compared to 50- to 75-µm isotropic resolutions reported in previous dMRI studies (15, 16, 58). Considering that the soma of mouse radial glial cells are ∼10 to 20 µm in diameter (60) and their processes can extend over 200 µm (60), our method offers an opportunity to depict mesoscopic cellular organization starting in midgestation with the radial organization of the telencephalic wall of E10.5 and E11.5 embryos (Fig. 1) that was not available in previous dMRI studies. Moreover, the 3D dMRM data with isotropic resolution also allows tracing of the fiber tracts from the very beginning at E12.5 and following their development during gestation. These rich dMRM contrasts could be potentially combined with new MRI contrasts (7, 61) and optical-imaging techniques (62, 63) to build a more-encompassing toolkit for studying cell migration, differentiation, and axon tract formation in the developing brain.

Developmental Continuum Reconstructed Using Computational Anatomy Techniques.

Techniques used here to characterize the dynamic processes occurring in brain development, including diffeomorphic mapping and quantification of volumetric deformation, have been developed under the framework of computational anatomy (19, 20). Applying these techniques to MRI data of developing human and animal brains have revealed whole-brain growth patterns and contributed to defining the growth trajectories of major brain structures (27, 28, 64, 65). The reconstruction of the 4D GVF by concatenating diffeomorphic mappings among consecutive developmental stages from E11.5 to E15.5 enables us to visualize the spatiotemporal changes of the embryonic mouse brain assembled in a common coordinate system.

Searching for structures within the 4D continuum that share similar developmental trajectories, we used a data-driven approach and identified five developmental modes with distinct trajectories. A similar approach has been recently used to produce adult mouse brain parcellation based on transcriptomic data (66). To understand the genetic driving force of these modes, we further performed a joint analysis with the rich genetic database in ADMBA utilizing the AGEA tool. Note that due to developmental difference () between the CD1 embryos used here and the C57BL/6J embryos used in ADMBA (with gene-expression data in one hemisphere and slide slices missing), we did not perform coregistration between the Allen ISH images and dMRM.

Four-Dimensional Continuum Connects Upstream Genetic Events with Downstream Neuroanatomy.

A major advantage of the proposed 4D continuum lies in its ability to colocalize and transform the genetic expression or neuroanatomical data along the temporal axis in pursuit of identifying their interactions not only at the concurrent but at distant time points as well. This is useful, because the observed anatomical phenotype may be a consequence of earlier or accumulated patterns of gene expression and related cellular events. We illustrated this concept with two examples: 1) mapping Slit2/Robo1 expressions at E12.5 to E15.5 when the emergence of commissural fibers was observed in dMRM (Fig. 5), and 2) tracing the origin of the cortical abnormality detected at E15.5 back to E12.5, when abnormal migration occurred in the Sas-4 model (Fig. 6). Note that the transformation between different stages assumed that brain topology remained relatively consistent across stages, which is generally valid between E11.5 and E15.5 (1) ().

However, the gene–neuroanatomy correlation analyzed in the current study only analyzed a small number of genes as a proof of concept, which limited the use of advanced bioinformatics tools (67–69) to explore large-scale molecular profiles for C57BL/6 or CD1 embryos (70, 71). In future studies, the dMRM-based developmental continuum could be integrated with spatially resolved transcriptomics (72) to link single-cell transcriptomic profiles with neuroanatomy in normal developing and genetically modified mouse brains.

In summary, we established a developmental continuum of the embryonic mouse brain based on ultra high-resolution dMRM and characterized the 4D neuroanatomy using integrated computational analysis, which revealed distinct modes of development in E10.5 through E15.5 embryos. By investigating the spatiotemporal changes of gene-expression patterns in the developing continuum, this tool can be used to investigate mechanisms controlling the development of mesoscopic neuroanatomy and opens the possibility of systematically studying the genetic control of brain development.

Materials and Methods

Embryonic Mouse Brain Specimens.

All animal procedures were approved by the Animal Use and Care Committee at the Johns Hopkins School of Medicine. Embryos were obtained by euthanizing time-mated, pregnant CD1 mice (Charles River Laboratories, Inc.). A total of 30 embryos at six gestational stages at embryonic day 10.5 (E10.5), E11.5, E12.5, E13.5, E14.5, and E15.5 (n = 5 per stage) were collected. The whole body of the E10.5 embryos and the brains of the E11.5 through E15.5 embryos were immersion fixed in 4% paraformaldehyde for over 1 wk and then transferred to phosphate buffered saline with 1 mM Gd-DTPA (Berlex Imaging, Wayne, NJ) for over 72 h to enhance the MR signals. During imaging, the specimens were embedded in 5% agarose gel (Sigma-Aldrich, St. Louis, MO) and placed in 10-mm diameter NMR tubes filled with Fomblin (Solvay Solexis, Thorofare, NJ) to prevent dehydration and susceptibility artifacts.

Gene Expression and ISH.

The genetic probes and ISH procedures used in this study have been previously described in ref. 73.

Diffusion MRI Acquisition.

All imaging was performed on a vertical 11.7-Tesla scanner (Bruker Biospin, Billerica, MA) with a birdcage volume coil (10-mm inner diameter), a Micro5 gradient system (maximum gradient strength of 3,000 mT/m), and a temperature control system. The E10.5 embryos were scanned using a solenoid coil, and E11.5 through E15.5 mouse brains were scanned using the 10-mm volume coil.

Three-dimensional high-resolution diffusion MRI was acquired using an in-house 3D diffusion-weighted gradient and spin-echo (DW-GRASE) sequence (35) with double-sampled echo-planar readout and twin navigator echoes for motion and phase corrections. The E10.5 embryos, E11.5 through E14.5 embryo brains, and E15.5 embryo brains were imaged at 30 μm, 35 μm, and 45 μm isotropic resolutions, respectively, with the following parameters: echo time = 30 ms, repetition time = 700 ms, two signal averages, diffusion gradient duration = 4 ms, diffusion gradient separation = 12 ms, 30 diffusion directions with a b-value of 1.2 ms/µm2 for E10.5, 1.3 ms/µm2 for E11.5 through E14.5, and 2.1 ms/µm2 for E15.5 brains. High b-values were used for later embryonic stages, as the brain diffusivity decreases with development (74). Note that due to the relaxation-induced MR signal decay during signal readout, the actual image resolution, as characterized by the point spread function (PSF), slightly differs from the nominal resolution set in the protocols. For example, the average full-width half-maximum of PSF for the E12.5 brain was 35.2 µm and 24.6 µm along the y and z phase-encoding directions, respectively ().

Diffusion Tensor and Fiber Orientation Density Reconstruction.

The 3D DW-GRASE data were reconstructed from k-space to images and zero padded to twice the raw image resolution in each dimension in MATLAB (MathWorks, Natick). Diffusion tensor (75) was reconstructed using the log-linear fitting method in DTI (Diffusion Tensor Imaging) Studio (http://www.mristudio.org), and the tensor-related metrics were generated, including the MD and FA. Three-dimensional reconstruction of white matter tracts was performed using the fiber assignment by continuous tracking method (76). Tracking was performed using manually defined regions of interest (ROIs) as seed regions with a FA threshold of 0.05 and a maximum angle of 60°. In order to resolve more-complex microstructural organizations, such as crossing fibers, we used the constraint spherical deconvolution method (77) to estimate the fiber orientation distributions in MRtrix (http://www.mrtrix.org). FOD-based, probabilistic, streamline tracking was performed to generate super-high–resolution TDIs at a grid-size 10-µm isotropic resolution.

Embryonic Mouse Brain Atlas Generation Based on dMRM.

Five embryonic brains at each embryonic stage were coregistered to generate population-averaged atlases per stage from E11.5 to E15.5 using an iterative procedure. Briefly, in the initial iteration, one representative embryonic brain was selected as the reference, and the remaining brains were registered to the reference brain, first by landmark-based rigid registration and then intensity-based affine registration and finally by nonlinear registration using a two-channel large deformation diffeomorphic metric mapping (LDDMM) (78) based on the mean diffusion-weighted image (mDWI) and FA contrasts, to obtain an averaged atlas brain. In the next iteration, all five brains were transformed to the atlas from the last iteration with affine and LDDMM registrations to obtain an updated atlas. The procedure was repeated three times to obtain the final atlas brain. The landmark, affine, and LDDMM registrations were performed using DiffeoMap (http://www.mristudio.org). The registration accuracy was accessed by independent sets of landmarks that were manually placed at both the brain contours and the white matter tracts on dMRM images (), and the distances between paired landmarks on the individual brains and atlases were 0.010 to 0.017 mm for E11.5 through E15.5, which was less than a voxel ().

Anatomical Development Quantified by Computational Anatomy.

In order to quantify the embryonic brain deformation with development, we performed longitudinal registration between consecutive embryonic stages. For instance, the E11.5 atlas was transformed to the E12.5 atlas, first by landmark-based rigid transformation and intensity-based rigid transformation and then by the two-channel LDDMM based on the mDWI and FA contrasts (79). In the cost function used by LDDMM, the term that controls the smoothness of the vector field was adjusted to capture macroscopic changes (>0.1 mm) in brain morphology. The cross-stage registration accuracy was accessed by independent sets of landmarks placed at the brain contours and the white matter tracts on atlas images (), and the landmark distances between the atlases were 0.24 to 0.42 mm for E11.5 through E15.5, which translated to one to two voxels ().

The voxelwise GVF between the adjacent stages was obtained from LDDMM-based transformation = [T], from which the Jacobian matrix () was calculated.

The strain tensor was obtained from as to represent the shape and primary orientation of the tensor. The Jacobian determinant JDet was also obtained from the and log transformed to quantify the extent of deformation, which characterizes the amount of brain growth from one embryonic day to the next. The dMRM atlas generated at intermediate stages (at 0.1-d interval from E11.5 to E15.5; and Movie S1) was obtained interpolating the GVF between adjacent stages.

Developmental Mode.

In order to analyze the dMRI metrics in a common space and obtain the developmental trajectories in a voxelwise manner, the MRI data at E11.5 through E15.5 were all transformed to the E13.5 atlas space by concatenating forward and inverse transformation matrices from the longitudinal registration as mentioned above (Fig. 2). Note that the E10.5 atlas was excluded in this analysis due to the apparent difference in topology compared with the other developmental stages. In the E13.5 common space, FA and J values were extracted from E11.5 to E15.5 as a feature vector [FAE11.5, FAE12.5, FAE13.5, FAE14.5, J,E11.5, J,E12.5, J,E13.5, J,E14.5] in a voxelwise matter. k-means clustering was performed on the feature vectors that clustered the voxels with similar developmental trajectories together. We varied the number of clusters from four to eight (), and a cluster number of six resulted in the most anatomically interpretable modes.

Registration between dMRM and ISH Images of Gene Expression.

Stacks of ISH slides of gene expression were registered to the gestational day–matched dMRM atlases with the following procedure. First, a small number of landmarks were placed on the ISH slides (seven to 11 per slide) and the corresponding mDWI slices for landmark-based affine registration to coarsely place the ISH data in the atlas space. Then, a large number of landmarks (200 to 300 per sample; ) were placed across the brain for landmark-based LDDMM transformation in DiffeoMap (http://www.mristudio.org). After this, the ISH colormaps were converted to gray-scale images and normalized for intensity-based registration with the mDWIs using LDDMM with mutual information as the cost function. The registration accuracy was accessed by independent sets of landmarks manually placed on the ISH images and dMRI atlas (), and the averaged distances between the paired landmarks were 0.043 ± 0.018 mm, 0.034 ± 0.025 mm, 0.027 ± 0.028 mm, and 0.025 ± 0.011 mm for E11.5, E12.5, E14.5, and E15.5, respectively, indicating that registration errors were within one to two voxels ().

Two-Dimensional Dendrogram of Gene Expression in Individual Developmental Modes.

The intensity of gene expression in the stained areas was extracted from the normalized ISH images within each developmental mode. In each mode, the log-transformed sum of expression intensities for the 10 coregistered genes at E11.5, E12.5, E14.5, and E15.5 formed a 2D matrix, and the intensities were standardized into [−1, 1]. Hierarchical clustering of the 2D matrices (80) was performed with unweighted, average linkage-clustering and Euclidean distance as the dissimilarity matrix in MATLAB.

Correlation Between Gene Expression and dMRM Metrics.

Based on the ADMBA parcellations, the log-transformed sum of gene-expression intensities in the prosomeric segments at E11.5, E12.5, E14.5, and E15.5 were obtained to form a 2D matrix for each gene, which was then normalized for comparison. Similarly, we calculated the averaged FA and mean/axial/radial diffusivities (MD/AD/RD) in each segment across gestation to obtain 2D spatiotemporal patterns. For each gene, we performed the 2D correlations between the expression intensity and the dMRM feature matrices to obtain the correlation coefficient and P values. Note that only 17 prosomeric segments in the telencephalon and diencephalon were used in this analysis as some of the ISH images did not cover the rostral brain regions.

dMRM of the Sas4 Embryonic Mouse Brains.

Ex vivo dMRM of wild-type and genetically modified Sas4 mouse brains was performed at embryonic day E15.5 with the same experimental setup. Sas-4 deletion led to a progressive loss of centrioles, resulting in detached RGPs from ventricular zone, and simultaneous removal of p53 could rescue RGP death and microcephaly but could not prevent RGP delocalization and abnormal neuronal migration (53). The Sas4 conditional allele was described in detail in Bazzi et al. (83), and the immunohistochemistry and confocal imaging procedures were described in detail in Insolera et al. (53).