We compute the full O(alpha_s) supersymmetric QCD corrections for neutralino-stop co-annihilation into electroweak gauge and Higgs bosons in the Minimal Supersymmetric Standard Model (MSSM). We show that these annihilation channels are phenomenologically relevant within the so-called phenomenological MSSM, in particular in the light of the observation of a Higgs-like particle with a mass of about 126 GeV at the LHC. We present in detail our calculation, including the renormalization scheme, the infrared treatment, and the kinematical subtleties to be addressed. Numerical results for the co-annihilation cross sections and the predicted neutralino relic density are presented. We demonstrate that the impact of including the corrections on the cosmologically preferred region of parameter space is larger than the current experimental uncertainty from WMAP data.