In this paper, we perform a full next-to-leading order (NLO) QCD calculation of neutralino scattering on protons or neutrons in the MSSM. We match the results of the NLO QCD calculation to the scalar and axial-vector operators in the effective field theory approach. These govern the spin-independent and spin-dependent detection rates, respectively. The calculations have been performed for general bino, wino and higgsino decompositions of neutralino dark matter and required a novel tensor reduction method of loop integrals with vanishing relative velocities and Gram determinants. Numerically, the NLO QCD effects are shown to be of at least of similar size and sometimes larger than the currently estimated nuclear uncertainties. We also demonstrate the interplay of the direct detection rate with the relic density when consistently analyzed with the program $\mathtt{DM@NLO}$.

For particle physics observables at colliders such as the LHC at CERN, it has been common practice for many decades 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. In this paper, we present the first quantitative study of the theoretical uncertainty of the neutralino dark matter relic density from scheme and scale variations. We first explain in detail how the renormalization scale enters the tree-level calculations through coupling constants, masses and mixing angles. We then demonstrate a reduction of the renormalization scale dependence through one-loop SUSY-QCD corrections in many different dark matter annihilation channels and enhanced perturbative stability of a mixed on-shell/$\bar{\rm DR}$ renormalization scheme over a pure $\bar{\rm DR}$ scheme in the top-quark sector. In the stop-stop annihilation channel, the Sommerfeld enhancement and its scale dependence are shown to be of particular importance. Finally, the impact of our higher-order SUSY-QCD corrections and their scale uncertainties are studied in three typical scenarios of the phenomenological Minimal Supersymmetric Standard Model with eleven parameters (pMSSM-11). We find that the theoretical uncertainty is reduced in many cases and can become comparable to the size of the experimental one in some scenarios.