The quadrupole collective Hamiltonian, based on relativistic energy density functionals, is extended to include a pairing collective coordinate. In addition to quadrupole shape vibrations and rotations, the model describes pairing vibrations and the coupling between shape and pairing degrees of freedom. The parameters of the collective Hamiltonian are determined by constrained self-consistent relativistic mean-field plus Bardeen-Cooper-Schrieffer (RMF+BCS) calculations in the space of intrinsic shape and pairing deformations. The effect of coupling between shape and pairing degrees of freedom is analyzed in a study of low-energy spectra and transition rates of four axially symmetric $N=92$ rare-earth isotones. When compared to results obtained with the standard quadrupole collective Hamiltonian, the inclusion of dynamical pairing increases the moment of inertia, lowers the energies of excited $0^+$ states and reduces the E0-transition strengths, in better agreement with data.