One of the most important mechanisms at the Large Hadron Collider (LHC) for the production of the pseudoscalar Higgs boson of the Minimal Supersymmetric Standard Model (MSSM) is the loop-induced gluon fusion process gg → A. The higher-order QCD corrections have been obtained a long time ago and turned out to be large. However, the genuine supersymmetric (SUSY–)QCD corrections have been obtained only in the limit of large SUSY particle masses so far. We describe our calculation of the next-to-leading-order (NLO) SUSY-QCD results with full mass dependence and present numerical results for a few representative benchmark points. We also address the treatment of the effective top and bottom Yukawa couplings, in the case of heavy SUSY particles, in terms of effective low-energy theories where the heavy degrees of freedom have been decoupled. Furthermore, we include a discussion of the relation between the SUSY-QCD corrections that we have computed and the Adler-Bardeen theorem for the axial anomaly. In addition, we apply our results to the gluonic and photonic pseudoscalar Higgs decays A → gg, γγ at NLO.