  We present the full $\mathcal{O}(\alpha_s)$ supersymmetric QCD corrections
for gaugino annihilation and co-annihilation into light and heavy quarks in the
Minimal Supersymmetric Standard Model (MSSM). We demonstrate that these
channels are phenomenologically relevant within the so-called phenomenological
MSSM. We discuss selected technical details such as the dipole subtraction
method in the case of light quarks and the treatment of the bottom quark mass
and Yukawa coupling. Numerical results for the (co-)annihilation cross sections
and the predicted neutralino relic density are presented. We show that the
impact of including the radiative corrections on the cosmologically preferred
region of the parameter space is larger than the current experimental
uncertainty from Planck data.
