The Laplace TF of a ZOH is usually taken as \$\small \dfrac{1-e^{-sT}}{s}\$, where \$\small T\$ is the sampling increment. Multiply this by \$\small G_p(s)\$, giving $$\small G^*_p(s)=(1-e^{-sT}) \frac{G_p(s)}{s}$$
Since \$\small e^{-sT}\$ and \$\small z^{-1}\$ represent pure time delays of \$\small T\$ in the s and z domains, respectively, we may write:
$$ \small 1-e^{-sT} \rightarrow \dfrac{z-1}{z}$$
The z-transform of the remaining term, \$\small \dfrac{G_p(s)}{s}\$, may be obtained from the tables {http://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html}
Working (very) approximately from your equations to check the result, the original Laplace TF may be simplified to: $$\small G_p(s)\approx\dfrac{ \omega_n^2}{s^2+\omega_n^2}$$
with \$\small \zeta=0.017\approx0\$; \$\small \omega_n=7\times 10^3rad/s\$; \$\small \omega_n\:T=0.547rad\$ \$\small\approx30^0\$.
Hence, after z-transforming \$\small \dfrac{G_p(s)}{s}\$ and multiplying by \$\small \dfrac{z-1}{z}\$, the approximate z-TF from the tables is: $$\small G^*_p(z)\approx\frac{0.13z+0.13}{z^2-1.73z+1}$$
Which compares favourably with the answer you give.