  For a given bounded domain $\Omega$ with smooth boundary in a smooth
Riemannian manifold $(\mathcal{M},g)$, we establish a procedure to get all the
coefficients of the asymptotic expansion of the trace of the heat kernel
associated with the perturbed polyharmonic Dirichlet-to-Neumann operator
$\Lambda_m$ ($m\ge 1$) as $t\to 0^+$. We also explicitly calculate the first
four coefficients of this asymptotic expansion. These coefficients (i.e., heat
invariants) provide precise information for the area and curvatures of the
boundary $\partial \Omega$ in terms of the spectrum of the perturbed
polyharmonic Steklov problem. In particular, when $m=1$ and $q\equiv 0$ our
work recovers the previous corresponding results in \cite{PS} and \cite{Liu3}.
