  We report a systematic study of nuclear matrix elements (NMEs) in
neutrinoless double-beta decays with a state-of-the-art beyond mean-field
covariant density functional theory. The dynamic effects of particle-number and
angular-momentum conservations as well as quadrupole shape fluctuations are
taken into account with projections and generator coordinate method for both
initial and final nuclei. The full relativistic transition operator is adopted
to calculate the NMEs. The present systematic studies show that in most of the
cases there is a much better agreement with the previous non-relativistic
calculation based on the Gogny force than in the case of the nucleus $^{150}$Nd
found in Song et al. [Phys. Rev. C 90, 054309 (2014)]. In particular, we find
that the total NMEs can be well approximated by the pure axial-vector coupling
term with a considerable reduction of the computational effort.
