How an equation in s domain is converted into z domian using zero hold method
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 is: $$\small G^*_p(z)\approx\frac{0.13z+0.13}{z^2-1.73z+1}$$

Which compares favourably with the answer you give.