Nearly 20 years ago Masatoshi Nei (1987) asserted that the study of evolution is focused on two types of questions: reconstructing the histories of individual groups of organisms and understanding the mechanisms that drive evolutionary change. Nei's dichotomy, which he used to introduce his book on Molecular Evolutionary Genetics, served as a rhetorical device around which he integrated emerging knowledge of genetic diversity at the molecular level, concluding that the analysis of DNA and protein data would rapidly extend knowledge on both fronts. But more important for Nei was the possibility that phylogenetic and evolutionary inference might be fused through the quantitative analysis of molecular data, and his book was a bold attempt to achieve that synthesis. Now, two decades later, the power of Nei's molecular genetic synthesis can hardly be doubted. The expansion of molecular studies has far outstripped any other evolutionary genre, including palaeontology. Population-level data on geographic variation in gene sequences, for example, provide insights into mutation rates, demography, migration, speciation, phylogenetic diversification, times of common ancestry, and response to geologic-scale environmental restructuring (Patton and Smith 1989; Avise 1994, 2000; Wakeley and Hey 1998; Hewitt 1999, 2000; Knowles 2004). The simultaneous inference of both pattern and process from genetic data is now the norm.

Nei's dichotomy still occurs in palaeontology, however, where quantitative synthesis of evolutionary process and phylogenetics is still in its infancy. While the field is rich, both with studies of the histories of groups of organisms and studies of evolutionary mechanisms, the data and methods employed in the two types of study are radically different. Discrete-character, parsimony-based cladistics remains the dominant paradigm for phylogeny reconstruction in palaeontology (Adrain 2001; Padian et al. 1994), but quantitative analysis of morphological size and shape (Raup 1967; MacLeod and Rose 1993; Foote 1994, 1997, 1999; Jernvall et al. 1996; Smith 1998; Gingerich and Uhen 1998; Smith and Lieberman 1999; Roopnarine 2001) or of taxonomic abundance and extinction (Raup 1976; Raup and Marshall 1980; Sepkoski 1997; Alroy 2003; Janis et al. 2004) are the most common palaeontological approaches to the study of evolutionary mechanisms. Indeed, many palaeontologists share a widespread mistrust that quantitative measures of morphology do not carry recoverable information about phylogenetic and evolutionary history (see reviews in MacLeod 2002; Humphries 2002; Felsenstein 2002).

One barrier to the quantitative synthesis of evolutionary mechanisms and phylogenetics in palaeontology is the lack of tools for comparing the evolutionary behaviour of complex morphological traits (e.g., skulls, teeth, lophophores, etc.) on both micro- and macroevolutionary timescales. How do we expect variation in morphology, which is shaped by development, growth, or functional constraints, etc., to scale between populations and phylogenies? Can we reconstruct population-level processes from comparative data drawn from a variety of related taxa or from a long stratigraphic sequence? How do we expect inter-taxon patterns of diversity in complex morphologies to be affected by morphological integration, or by different modes of selection, such as directional, stabilizing, or drift? Progress has been made in regards to some of these issues for simple morphological measures, such as size and ratios (Gingerich 1993, 2001; Wagner 1996; Roopnarine 2001, 2003; Polly 2001a). The simulation method presented here should help make further progress towards a morphological synthesis.

The simulations presented here make it clear that certain aspects of microevolutionary process can be detected by comparison of morphologies separated by long stretches of phylogenetic time. The scaling between multivariate morphological distance and the duration of phylogenetic divergence is linear under directional selection, curvilinear under randomly varying selection or drift, and stochastically constant under stabilizing selection. Comparison of the teeth of different populations, species, and genera of shrew demonstrates that the observed scaling between distance and divergence is similar to that under fluctuating selection.

Furthermore, comparison of empirical data to simulation results demonstrates that molar tooth form does not evolve by purely genetic drift. Drift is a function of morphological variance and population size, and is more significant at very small effective population sizes. Even when the mean effective population size is assumed to be only 70 individuals, drift does not produce enough change over even millions of generations to explain the observed morphological divergences in shrew molar shape. Consequently, the mode of evolution in shrew molar morphology is likely selection that fluctuates from generation to generation in intensity and direction.

The simulations also indicate that morphometric shape is likely to have a strong phylogenetic component that can be utilized for phylogeny reconstruction. When the simulation is given phylogenetic structure, related tip taxa share visible derived similarity. Closely related taxa are also located close to one another in principal components space, locations that can be estimated using principal components analysis of the tip taxa. Maximum-likelihood methods for continuous traits (Felsenstein 1973) are appropriate for estimating trees from such data because they are based on the assumption of Brownian-motion evolution, which is equivalent to the randomly fluctuation selection patterns found here. Previous empirical studies of morphometric-based phylogenetic trees have found that maximum-likelihood trees based on shape data correspond well to the patterns of relationship suggested by independent data (Polly 2003a, 2003b).

The simulation method presented here could be elaborated to incorporate more complicated evolutionary assumptions. The simulations run in this paper were based on P, the phenotypic variance-covariance matrix; however, the same method could be applied to G, the additive genetic covariance matrix, if it is available. Furthermore, the simulations could be run using more sophisticated quantitative genetic models. The simplest of these improvements would be to allow P to evolve, but more complicated covariance models that incorporate developmental interactions, as well as genetic models could also be used (Johnson and Porter 2001; Wolf et al. 2001; Rice 2002; Salazar-Cuidad and Jernvall 2002).