© 2021 by the authors. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.An adjoint-based aerodynamic shape optimization framework for natural laminar flows is developed. The laminar to turbulent-transition onset is predicted by the correlation-based Bas–Cakmakcioglu transition model that is coupled with the Spalart–Allmaras turbulence model. A discrete adjoint implementation is subsequently developed. Automatic differentiation is utilized to construct the partial derivatives in the discrete adjoint formulation. The turbulence and transition models are fully coupled into the sensitivity derivative evaluations, as well as into the objective function evaluations. The sensitivity derivatives evaluated by the discrete adjoint solver are validated against those of the finite central difference method. The discrete adjoint-based aerodynamic shape optimization framework developed for natural laminar flows is successfully employed to optimize aerodynamic characteristics of the NLF(1)-0416 airfoil and a low-aspect-ratio wing.