  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}$.
