We describe production of heavy quarkonia in pA collisions within the dipole approach, assuming dominance of the perturbative color-singlet mechanism (CSM) in the pT -integrated cross section. Although accounting for a nonzero heavy Q-Q separation is a higher twist correction, usually neglected, we found it to be the dominant source of nuclear effects, significantly exceeding the effects of leading twist gluon shadowing and energy loss. Moreover, this contribution turns out to be the most reliably predicted, relying on the precise measurements of the dipole cross section at HERA. The nuclear suppression of quarkonia has been anticipated to become stronger with energy, because the dipole cross section steeply rises. However, the measured nuclear effects remain essentially unchanged within the energy range from RHIC to the LHC. A novel production mechanism is proposed, which enhances the charmonium yield. Nuclear effects for the production of J/ψ, ψ(2S), Υ(1S) and Υ(2S) are calculated, in agreement with data from RHIC and LHC. The dipole description offers a unique explanation for the observed significant nuclear suppression of ψ(2S) to J/ψ ratio, related to the nontrivial features of the ψ(2S) wave function.