To gain for a more precise prediction of the supersymmetric dark matter relic density, we have calculated next-to-leading order SUSY-QCD corrections to neutralino (co)annihilation processes including Coulomb enhancement effects. We demonstrate that these corrections can have significant impact and are thus of general interest for parameter studies. For particle physics observables at colliders, it is common practice to estimate the theoretical uncertainty by studying the variations of the predicted cross sections with a priori unpredictable scales. In astroparticle physics, this has so far not been possible, since most of the observables were calculated at Born level only, so that the renormalization scheme and scale dependence could not be studied in a meaningful way. Thus, we present the fi rat quantitative study of the theoretical uncertainty of the neutralino dark matter relic density from scheme and scale variations.