Background: The coexistence of different geometric shapes at low energies presents a universal structure phenomenon that occurs over the entire chart of nuclides. Studies of the shape coexistence are important for understanding the microscopic origin of collectivity and modifications of shell structure in exotic nuclei far from stability. Purpose: The aim of this work is to provide a systematic analysis of characteristic signatures of coexisting nuclear shapes in different mass regions, using a global self-consistent theoretical method based on universal energy density functionals and the quadrupole collective model. Method: The low-energy excitation spectrum and quadrupole shape invariants of the two lowest 0+ states of even-even nuclei are obtained as solutions of a five-dimensional collective Hamiltonian (5DCH) model, with parameters determined by constrained self-consistent mean- field calculations based on the relativistic energy density functional PC-PK1, and a finite-range pairing interaction. Results: The theoretical excitation energies of the states, 2+1, 4+1, 0+2, 2+2, 2+3, as well as the B(E2 ; 0+1→2+1) values, are in very good agreement with the corresponding experimental values for 621 even-even nuclei. Quadrupole shape invariants have been implemented to investigate shape coexistence, and the distribution of possible shape-coexisting nuclei is consistent with results obtained in recent theoretical studies and available data. Conclusions: The present analysis has shown that, when based on a universal and consistent microscopic framework of nuclear density functionals, shape invariants provide distinct indicators and reliable predictions for the occurrence of low-energy coexisting shapes. This method is particularly useful for studies of shape coexistence in regions far from stability where few data are available.