This work proposes a framework for fully-automatic gradient-based constrained aerodynamic shape optimization in a multistage turbomachinery environment. A turbomachinery solver which solves the Reynolds-averaged Navier-Stokes (RANS) equations to a steady-state in both rotating and stationary domains is developed. Characteristic-based inlet and outlet boundary conditions are imposed, while adjacent rotor and stator rows are coupled by mixing-plane interfaces. To allow for an efficient but accurate gradient calculation, the turbomachinery RANS solver is adjointed at a discrete level. The systematic approach for the development of the discrete adjoint solver is discussed. Special emphasis is put on the development of the turbomachinery specific features of the adjoint solver, i.e. on the derivation of flow-consistent adjoint inlet and outlet boundary conditions and, to allow for a concurrent rotor-stator optimization and stage coupling, on the development of an exact adjoint counterpart to the non-reflective, conservative mixing-plane formulation used in the flow solver.The adjoint solver is validated by comparing its sensitivities with finite-difference gradients obtained from the flow solver. A parallelized, automatic grid perturbation scheme utilizing radial basis functions, which is accurate and robust as well as able to handle complex multi-block grid configurations, is employed to calculate the gradient