Accelerating the evaluation of analytical nuclear energy gradients in wavefunction-based quantum chemical methods and extending these approaches to large molecular systems remain challenging tasks. In this work, we developed analytical nuclear gradients for second-order Møller-Plesset perturbation theory (MP2) employing pair-natural orbitals (PNOs) derived from orthonormal localized virtual molecular orbitals (LVMOs). In contrast to projected atomic orbital (PAO)-based formulations, our approach introduces two key distinctions. First, constraints arising from the LVMO construction are incorporated into the Lagrangian, which necessitates solving an additional Z-vector equation, namely the coupled-perturbed virtual localization (CP-VL) equation. Second, due to the orthonormal nature of the LVMOs, the derivatives of the PNO coefficients with respect to both MO coefficients and nuclear coordinates vanish, which simplifies the formulation of analytical energy gradients. Additionally, we developed an automatic derivation program to alleviate the technical complexity associated with deriving gradient expressions within the PNO framework. We also briefly report an extension of this framework to analytical nuclear gradients for the CC2 method using PNOs based on LVMOs for ground states.