There is a common belief that the main uncertainties in the theoretical analysis of neutrinoless double beta ($0\nu\beta\beta$) decay originate from the nuclear matrix elements. Here, we uncover another previously overlooked source of potentially large uncertainties stemming from non-perturbative QCD effects. Recently perturbative QCD corrections have been calculated for all dimension 6 and 9 effective operators describing $0\nu\beta\beta$-decay and their importance for a reliable treatment of $0\nu\beta\beta$-decay has been demonstrated. However, these perturbative results are valid at energy scales above $\sim 1$ GeV, while the typical $0\nu\beta\beta$-scale is about $\sim 100$ MeV. In view of this fact we examine the possibility of extrapolating the perturbative results towards sub-GeV non-perturbative scales on the basis of the QCD coupling constant "freezing" behavior using Background Perturbation Theory. Our analysis suggests that such an infrared extrapolation does modify the perturbative results for both short-range and long-range mechanisms of $0\nu\beta\beta$-decay in general only moderately. We also discuss that the tensor$\otimes$tensor effective operator can not appear alone in the low-energy limit of any renormalizable high-scale model and then demonstrate that all five linearly independent combinations of the scalar and tensor operators, that can appear in renormalizable models, are infrared stable.