The layer-by-layer deposition process used in material extrusion (ME) additive manufacturing results in inter- and intra-layer bonds that reduce the mechanical performance of printed parts. Multi-axis (MA) ME techniques have shown potential for mitigating this issue by enabling tailored deposition directions based on loading conditions in three dimensions (3D). Planning deposition paths leveraging this capability remains a challenge, as an intelligent method for assigning these directions does not exist. Existing literature has introduced topology optimization (TO) methods that assign material orientations to discrete regions of a part by simultaneously optimizing material distribution and orientation. These methods are insufficient for MA–ME, as the process offers additional freedom in varying material orientation that is not accounted for in the orientation parameterizations used in those methods. Additionally, optimizing orientation design spaces is challenging due to their non-convexity, and this issue is amplified with increased flexibility; the chosen orientation parameterization heavily impacts the algorithm’s performance. Therefore, the authors (i) present a TO method to simultaneously optimize material distribution and orientation with considerations for 3D material orientation variation and (ii) establish a suitable parameterization of the orientation design space. Three parameterizations are explored in this work: Euler angles, explicit quaternions, and natural quaternions. The parameterizations are compared using two benchmark minimum compliance problems, a 2.5D Messerschmitt–Bölkow–Blohm beam and a 3D Wheel, and a multi-loaded structure undergoing (i) pure tension and (ii) three-point bending. For the Wheel, the presented algorithm demonstrated a 38% improvement in compliance over an algorithm that only allowed planar orientation variation. Additionally, natural quaternions maintain the well-shaped design space of explicit quaternions without the need for unit length constraints, which lowers computational costs. Finally, the authors present a path toward integrating optimized geometries and material orientation fields resulting from the presented algorithm with MA–ME processes.