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