© 2018 authors. Published by the American Physical Society. There is a common belief that the main uncertainties in the theoretical analysis of neutrinoless double beta (0νββ) decay originate from the nuclear matrix elements. Here, we uncover another previously overlooked source of potentially large uncertainties stemming from nonperturbative QCD effects. Recently perturbative QCD corrections have been calculated for all dimension 6 and 9 effective operators describing 0νββ-decay and their importance for a reliable treatment of 0νββ-decay has been demonstrated. However, these perturbative results are valid at energy scales above ∼1 GeV, while the typical 0νββ scale is about ∼100 MeV. In view of this fact we examine the possibility of extrapolating the perturbative results towards sub-GeV nonperturbative 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νββ-decay in general only moderately. We also discuss that the tensor - tensor effective operator cannot 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, which can appear in renormalizable models, are infrared stable.