  We extend the multi-reference covariant density functional theory (MR-CDFT)
by including fluctuations in quadrupole deformations and average isovector
pairing gaps simultaneously for the nuclear matrix elements (NMEs) of
neutrinoless double-beta $(0\nu\beta\beta)$ decay in the candidate nuclei
$^{76}$Ge, $^{82}$Se, $^{100}$Mo, $^{130}$Te, and $^{136}$Xe assuming the
exchange of either light or heavy neutrinos. The results indicate a linear
correlation between the predicted NMEs and the isovector pairing strengths, as
well as the excitation energies of $2^{+}_1$ and $4^{+}_1$ states. By adjusting
the pairing strengths based on the excitation energies of the $2^{+}_1$ states,
we calculate the NMEs for $0\nu\beta\beta$ decay, which are reduced by
approximately $12\%$ to $62\%$ compared to the results obtained in the previous
studies by Song et al. [Phys. Rev. C95, 024305 (2017)]. Additionally, upon
introducing the average isovector pairing gap as an additional generator
coordinate in the calculation, the NMEs increase by a factor ranging from
$56\%$ to $218\%$.
