A hierarchical Bayesian framework accommodates intraspecific and interspecific variation in multivariate traits

Phylogenetic comparative methods are a critical tool in biology, providing the framework to test evolutionary hypotheses of phenotypic diversification. Accommodating intraspecific variation in these analyses is critical for accurate evolutionary inference, but current multivariate methods either assume traits evolve independently or that all taxa share the same intraspecific covariance structure. Violations of these assumptions can produce biased estimates of evolutionary parameters. Here, we introduce a hierarchical Bayesian framework for multivariate traits that jointly estimates taxon-specific intraspecific covariance structures alongside the underlying evolutionary process. This framework propagates uncertainty from sample size discrepancies and missing data, enabling the incorporation of highly variable morphological traits into phylogenetic analyses. Analyses of simulated data demonstrate that our framework achieves well-calibrated coverage (95%), whereas the standard practice of treating taxon means as known without error reduces coverage to 70%, confirming that ignoring intraspecific variation produces systematically biased evolutionary inference. We apply this framework to estimate the evolutionary rates and intraspecific distributions underlying perikymata distribution diversity in great apes, including modern humans and Neandertals. We show that, compared to other great apes, canine perikymata spacing in the genus Homo likely evolved under a substantially different regime than non-human apes, with cervical enamel evolving both rapidly and in a coordinated, modular fashion. We further find that Gorilla and Pongo show striking conservation in perikymata spacing relative to late Homo and Pan. These results and our method, which is applicable to other multivariate traits, provide the first phylogenetically rigorous characterization of an enamel growth trait across the great ape clade and establish intraspecific uncertainty propagation as a necessary component of multivariate phylogenetic analysis.