We argue that the inclusive $J/\psi$ production in proton-proton collisions might get significant contribution from fusion of three gluons. We analyze numerically the subprocesses with direct formation of $J/\psi$ from three gluons, and with formation of $J/\psi+g$ states. We found that numerically the three-gluon fusion gives a substantial contribution to the $J/\psi$ production, and is able to describe the experimentally observed shapes of the rapidity, momenta and multiplicity distributions. We also demonstrate that description of rapid multiplicity growth seen in experiment presents a challenge in approaches based on gluon-gluon fusion. For this reason we believe that the available multiplicity dependence data provide a strong evidence in favor of the 3-gluon mechanism of $J/\psi$ production.