The rare Higgs boson decay $H\to Z\gamma$ is forbidden at tree-level. In the Standard Model, it is loop-mediated through a $W$ boson or a heavy quark. We analytically compute the QCD correction to the heavy quark loop, confirming earlier purely numerical results, that were obtained for on-shell renormalization. The small quark mass expansion of the decay matrix element contains only single-logarithmic contributions at each perturbative order, which is in contrast to the double logarithms observed in $H\to \gamma\gamma$. We investigate the numerical interplay of bottom and top quark contributions, and the dependence of the result on the renormalization scheme.