Engineers and scientists are increasingly relying on high-fidelity numerical simulations. Within these simulations, mesh adaptation is useful for obtaining accurate predictions of an output of interest subject to a computational cost constraint. In the quest for accurately predicting outputs in problems with time-dependent solution features, a fully unstructured coupled spacetime approach has been shown to be useful in reducing the cost of the overall simulation. However, for the simulation of unsteady three-dimensional partial differential equations (PDEs), a four-dimensional mesh adaptation tool is needed. This work develops the first anisotropic metric-conforming four-dimensional mesh adaptation tool for performing adaptive numerical simulations of unsteady PDEs in three dimensions. The theory and implementation details behind our algorithm are first developed alongside an algorithm for constructing four-dimensional geometry representations. We then demonstrate our algorithm on three-dimensional benchmark cases and it appears to outperform existing implementations, both in metric-conformity and expected tetrahedra counts. We study the utility of the mesh adaptation components to justify the design of our algorithm. We then develop four-dimensional benchmark cases and demonstrate that metric-conformity and expected pentatope counts are also achieved. This is the first time anisotropic four-dimensional meshes have been presented in the literature. Next, the entire mesh adaptation framework, Mesh Optimization via Error Sampling and Synthesis (MOESS), is extended to the context of finding the optimal mesh to represent a function of four variables. The mesh size and aspect ratio distributions of the optimized meshes match the analytic ones, thus verifying our framework. Finally, we apply MOESS in conjunction with the mesh adaptation tool to perform the first four-dimensional anisotropic mesh adaptation for the solution of the advection-diffusion equation. The optimized meshes effectively refine the solution features corresponding to both a boundary layer solution as well as an expanding spherical wave.