The possibility of taking into account unsteady flow effects if performing turbomachinery shape optimization is attractive to accurately address inherently time dependent design problems. The harmonic balance method is an efficient solution for computational fluid dynamics problems of turbomachinery characterized by quasi-periodic flows. If applied in combination with adjoint methods, it enables the possibility to deal with unsteady fluid-dynamic design in a cost effective manner, opening the way towards multi-disciplinary applications. This paper presents the development of a novel fully-turbulent discrete adjoint based on the time domain harmonic balance method and its application to the constrained fluid dynamic optimization of an axial turbine stage. As opposed to previous works, the proposed method does not require any assumption on the turbulent eddy viscosity and on the set of input frequencies. The results show that the method provides accurate gradients, if compared with second order finite differences, and significant deviation with respect to the sensitivity computed with the constant eddy viscosity approximation. The application of the method to the fluid-dynamic shape optimization of the exemplary stage leads to improve the total-to-static efficiency of 0.8%. The efficiency increase is found to be higher than that obtained by means of a steady state optimization method.