A slight modification of Anton Petrunin's answer, namely Müller's idea brought up in Igor Belegradek's comment, allows one to deduce the desired result starting from the classical statement of the isometric embedding theorem, i.e.

Any Riemannian manifold $(M^n,g)$ admits a smooth isometric embedding into a Euclidean space.

Take again a smoothing $\phi$ of $\frac{1}{4}\operatorname{dist}_p$ with $|\phi-\frac{1}{4}\operatorname{dist}_p|\le 1$ and $|\nabla\phi|_g\le\frac{1}{2}$. This can be done by taking a locally finite open cover $\{U_j\}$ of $M$ such that $\overline{U}_j$ is compact and lies in a coordinate chart, together with a subordinated partition of unity $\{\rho_j\}$; since $|\nabla\operatorname{dist}_p|_g\le 1$ a.e., we can approximate $\frac{1}{4}\operatorname{dist}_p$ on $U_j$ by a smooth function $\phi_j$ with $$|\nabla\phi_j|_g\le\frac{3}{8},\qquad|\phi_j-\frac{1}{4}\operatorname{dist}_p|\le\min\bigg\{1,\frac{1}{8N_j'\|\nabla\rho_j\|_{L^\infty}}\bigg\},$$ where $N_k:=\#\{\ell:U_k\cap U_\ell\neq\emptyset\}$ and $N_j':=\max\{N_k\mid U_j\cap U_k\neq\emptyset\}$; finally, $\phi:=\sum\rho_j\phi_j$ has $|\phi-\frac{1}{4}\operatorname{dist}_p|\le 1$ and, being $\nabla\phi=\sum_k\rho_k\nabla\phi_k+\sum_k\nabla\rho_k(\phi_k-\frac{1}{4}\operatorname{dist}_p)$, on each $U_j$ we have $|\nabla\phi|_g\le\frac{3}{8}+\sum_{k:U_j\cap U_k\neq\emptyset}\frac{1}{8N_k'}\le\frac{1}{2}$ (as $N_k'\ge N_j$).

Now apply the classical Nash embedding theorem to the Riemannian manifold $(M^n,g-d\phi\otimes d\phi)$, obtaining an isometric embedding $F$. Then $(F,\phi)$ is the desired proper isometric embedding of $(M^n,g)$, since $\operatorname{dist}_p$ (and thus $\phi$) is a proper map by completeness.

A slight modification of Anton Petrunin's answer, namely Müller's idea brought up in Igor Belegradek's comment, allows one to deduce the desired result starting from the classical statement of the isometric embedding theorem, i.e.

Any Riemannian manifold $(M^n,g)$ admits a smooth isometric embedding into a Euclidean space.

Fix $p\in M$ and take again a smoothing $\phi$ of $\frac{1}{4}\operatorname{dist}_p$ with $|\phi-\frac{1}{4}\operatorname{dist}_p|\le 1$ and $|\nabla\phi|_g\le\frac{1}{2}$ (see below for the details). Now apply the classical Nash embedding theorem to the Riemannian manifold $(M^n,g-d\phi\otimes d\phi)$, obtaining an isometric embedding $F$. Then $(F,\phi)$ is the desired proper isometric embedding of $(M^n,g)$, since $\operatorname{dist}_p$ (and thus $\phi$) is a proper map by completeness.

Construction of $\boldsymbol{\phi}$. Take a locally finite open cover $\{U_j\}$ of $M$ such that $\overline{U}_j$ is compact and lies in a coordinate chart, together with a subordinated partition of unity $\{\rho_j\}$. Since $|\nabla\operatorname{dist}_p|_g\le 1$ a.e., we can approximate $\frac{1}{4}\operatorname{dist}_p$ on $U_j$ by a smooth function $\phi_j$ with $$|\nabla\phi_j|_g\le\frac{3}{8},\qquad|\phi_j-\frac{1}{4}\operatorname{dist}_p|\le\min\bigg\{1,\frac{1}{8N_j'\|\nabla\rho_j\|_{L^\infty}}\bigg\},$$ where $N_k:=\#\{\ell:U_k\cap U_\ell\neq\emptyset\}$ and $N_j':=\max\{N_k\mid U_j\cap U_k\neq\emptyset\}$. Finally, $\phi:=\sum\rho_j\phi_j$ has $|\phi-\frac{1}{4}\operatorname{dist}_p|\le 1$ and, being $$\nabla\phi=\sum_k\rho_k\nabla\phi_k+\sum_k\nabla\rho_k(\phi_k-\frac{1}{4}\operatorname{dist}_p),$$ on each $U_j$ (and thus on $M$) we have $|\nabla\phi|_g\le\frac{3}{8}+\sum_{k:U_j\cap U_k\neq\emptyset}\frac{1}{8N_k'}\le\frac{1}{2}$ (as $N_k'\ge N_j$).