We study exclusive quarkonium production in the dipole picture at next-to-leading order (NLO) accuracy, using the non-relativistic expansion for the quarkonium wavefunction. This process offers one of the best ways to obtain information about gluon distributions at small $x$, in ultraperipheral heavy ion collisions and in deep inelastic scattering. The quarkonium light cone wave functions needed in the dipole picture have typically been available only at tree level, either in phenomenological models or in the nonrelativistic limit. In this paper, we discuss the compatibility of the dipole approach and the non-relativistic expansion and compute NLO relativistic corrections to the quarkonium light-cone wave function in light-cone gauge. Using these corrections we recover results for the NLO decay width of quarkonium to $e^{+}e^{-}$ and we check that the non-relativistic expansion is consistent with ERBL evolution and with B-JIMWLK evolution of the target. The results presented here will allow computing the exclusive quarkonium production rate at NLO once the one loop photon wave function with massive quarks, currently under investigation, is known.