We present the two-loop QCD corrections to the amplitude of the Higgs production associated with a $Z$ boson via the bottom quark-antiquark annihilation channel with a non-vanishing bottom-quark Yukawa coupling. The computation is performed by projecting the D-dimensional scattering amplitude directly onto a set of Lorentz structures related to the linear polarisation states of the $Z$ boson. We cross-check the finite remainders through a computation based on conventional form factor decomposition. We show that for physical observables, an ultimate D-dimensional form factor decomposition of amplitudes is not necessary which has a huge potential to simplify a multiloop computation. We compute numerically the resulting cross sections under the soft-virtual approximation to NNLO and find it three orders of magnitude smaller than that of the s-channel.