We investigate the energy dependence of the photo-production cross-section of vector mesons J/Ψ and ϒ, measured by both HERA experiments H1 and ZEUS in electron-proton collisions and by LHC experiments ALICE, CMS and LHCb in ultra-peripheral proton-proton and ultraperipheral proton-lead collisions. Our study uses 2 particular fits of inclusive unintegrated gluon distribution, based on non-linear Balitsky-Kovchegov evolution (Kutak-Sapeta gluon; KS) and next-to-leading order Balitsky-Fadin-Kuraev-Lipatov evolution (Hentschinski-Sabio Vera-Salas gluon; HSS). We find that linear next-to-leading order BFKL evolution can only describe production at highest energies, if perturbative corrections are increased to unnaturally large values; rendering this corrections small, the growth with energy is too strong in the LHC region and the description of J/Ψ data fails. For the KS gluon we find that an accurate description of J/Ψ data is possible if non-linear corrections to low x QCD evolution are taken into account; without such correction a description of data fails. We interpret this observation as a clear signal for the presence of high gluon densities in low x the proton, characteristic for the onset of gluon saturation.