In addition to shape oscillations, low-energy excitation spectra of deformed nuclei are also influenced by pairing vibrations. The simultaneous description of these collective modes and their coupling has been a long-standing problem in nuclear structure theory. Here we address the problem in terms of self-consistent mean-field calculations of collective deformation energy surfaces, and the framework of the interacting boson approximation. In addition to quadrupole shape vibrations and rotations, the explicit coupling to pairing vibrations is taken into account by a boson-number non-conserving Hamiltonian, specified by a choice of a universal density functional and pairing interaction. An illustrative calculation for $^{128}$Xe and $^{130}$Xe shows the importance of dynamical pairing degrees of freedom, especially for structures built on low-energy $0^+$ excited states, in $gamma$-soft and triaxial nuclei.