An important ingredient for the calculation of Higgs boson properties in the infinite top quark mass limit is the knowledge of the effective coupling between the Higgs bosons and gluons, i.e. the Wilson coefficients CH and CHH for one and two Higgs bosons, respectively. In this work we calculate for the first time CHH to four loops in a direct, diagrammatic way, discussing in detail all issues arising due to the renormalization of operator products. Furthermore, we also calculate the Wilson coefficient CH for the coupling of a single Higgs boson to gluons as well as all four loop decoupling relations in QCD with general SU(Nc) colour factors. The latter are related to CH and CHH via low-energy theorems, which are used to obtain five-loop results for the Wilson coefficients.