During the last two decades, there is a continuous effort to design wings for aircraft with improved efficiency and reduced environmental footprint. Among others, there is a focus on the so-called Natural Laminar Flow (NLF) wings by extending the area over them where the flow is laminar, leading to reduced drag (and emissions) and, as a consequence, to environmental sustainability. The design of NLF wings requires computational tools that can accurately capture transitional phenomena and, also, incorporate this information into the optimization loop. In aerodynamic shape optimization, gradient-based algorithms usually rely on the adjoint method to compute gradients, since this is the most cost-effective way to do so. In this work, the continuous adjoint method for transitional flows is developed; in specific, the continuous adjoint to the γ-Re transition model coupled with the Spalart-Allmaras turbulence model is formulated and presented for the first time. The developed method is used for the shape optimization of a High Aspect-Ratio Wing Business Jet for minimum drag coefficient under lift and thickness constraints. The in-house GPU-accelerated PUMA code is used to carry out the optimization. The adjoint equations and their boundary conditions will be presented; consistent discretization schemes for the adjoint equations will be discussed in a companion paper. This work is part of the NEXTAIR project funded by the European Union under Grant Agreement No. 101056732.