We provide exact solutions of axisymmetric neutral vortices, that is, vortices with zero amount of potential vorticity anomaly (PVA), in three-dimensional (3D) quasi-geostrophic (QG) flows with distributed vorticity and density stratification. These solutions are linear combinations of PVA spherical layer-modes, which consist of conveniently normalized spherical Bessel functions of order 0, truncated by a zero of the spherical Bessel function of order 1. It is shown that, depending on the superposition of the different spherical layer-modes, some vortices remain axisymmetrically robust to small amplitude PVA perturbations, while other vortices evolve to stable QG tripoles. The robust axisymmetric 3D QG vortices analyzed here do not generate exterior potential flow and may provide an explanation to the persistence of baroclinic eddies in the ocean.