PDE-constrained shape optimization problems are characterized by target functionals that depend both on the shape of a domain (the control) and on the solution of a boundary value problem formulated on that domain (the state). Such optimization problems are commonly solved with mesh-moving methods. First, an initial mesh is constructed. Then, a (first order) finite element approximation of the state is computed and used to evaluate the misfit functional and its shape derivative. Finally, the coordinates of the mesh nodes are update. This procedure is repeated until convergence. We describe a higher-order mesh-moving method based on deformation diffeomorphisms. This approach supports naturally the use of higher-order finite element approximations of the state.