In the virtual presence of a heavy quark t, the interactions of a CP-odd scalar boson A, with mass M A 2M t , with gluons and light quarks can be described by an effective Lagrangian. We analytically derive the coefficient functions of the respective physical operators to three loops in quantum chromodynamics (QCD), adopting the modified minimal-subtraction (MS) scheme of dimensional regularization. Special attention is paid to the proper treatment of the γ 5 matrix and the Levi-Civita ε tensor in D dimensions. In the case of the effective ggA coupling, we find agreement with an all-order prediction based on a low-energy theorem in connection with the Adler-Bardeen non-renormalization theorem. This effective Lagrangian allows us to analytically evaluate the next-to-leading QCD correction to the A -> gg partial decay width by considering massless diagrams. For M A = 100 GeV, the resulting correction factor reads 1 + (22112)α ( 5 ) s (M A )/π + 165.9 (α ( 5 ) s (M A )/π) 2 1 + 0.68 + 0.23. We compare this result with predictions based on various scale-optimization methods.