For both fluid mechanics and practical engineering problems, it is important to study the effect of non-parallelism on the natural transition in boundary layers around underwater axisymmetric bodies. This paper develops the method of harmonic linearized Navier–Stokes (HLNS) equations for such boundary layers. The method can be used to study the evolution of small disturbances considering the non-parallelism of the basic flow. Based on the results of the HLNS equations, the transition positions of the boundary layers are predicted by the eN method. Because traditional linear stability theory (LST) neglects non-parallelism, the difference between the results using the HLNS method and those given by LST represents the effects of non-parallelism. Numerical calculations are performed for five classical forebody shapes, and the effect of non-parallelism is identified. (i) At each streamwise location, non-parallelism suppresses low-frequency disturbances while promoting high-frequency ones. The influence on high-frequency disturbances is more obvious. (ii) The effect of non-parallelism on the neutral curve focuses on the region near the critical instability position. Non-parallelism slightly delays the critical instability position for most forebody shapes and pushes the unstable zone toward the high-frequency direction while broadening its frequency bandwidth. (iii) The streamwise range at each growth-rate level is widened, implying that non-parallelism destabilizes the boundary layers. (iv) The transition occurs earlier, indicating that non-parallelism promotes the transition. (v) Non-parallelism shifts the dangerous frequency band toward the high-frequency direction. (vi) Non-parallelism enhances the wall pressure fluctuations (WPF) in the unstable laminar zone.