The use of high-fidelity computational fluid dynamics (CFD) tools in turbomachinery design has seen a continuous increase as a result of computational power growth and numerical methods improvement. These tools are often used in optimization environments, where gradient-based optimization algorithms are the most common due to their efficiency. In cases where the optimization contains a large number of design variables, the adjoint approach for calculating the gradients is beneficial, as it provides a way of obtaining function sensitivities with a computational cost that is nearly independent of the number of design variables. The interaction between adjacent blade rows is of utmost importance for the performance of multistage turbomachines. The most commonly used method to address these effects (i.e. coupling in the simulation of multiple rows) is the mixing-plane treatment, that has become a standard industrial tool in the design environment. In this paper, the formulation and implementation of an adjoint solver for multistage turbomachinery applications are presented, namely the adjoint counterpart of the mixing-plane formulation used in the direct solver. The solver is developed using the discrete ADjoint approach, where the partial derivatives required for the assembly of the adjoint system of equations are obtained using automatic differentiation tools. The sensitivity of several performance metrics relative to neighbor blade/hub row geometry and boundary conditions are shown to highlight the physical coupling in multi-row turbomachines. The verification of the adjoint multistage solver against the finite-difference approach is performed successfully with relative differences below 1 %.