We may use the Hankel transform representation of the Gaussian term (DLMF):
\begin{equation}
\int_{0}^{\infty}J_{\nu}\left(bt\right)\exp\left(-p^{2}t^{2}\right)t^{\nu+1}%
\mathrm{d}t=\frac{b^{\nu}}{(2p^{2})^{\nu+1}}\exp\left(-\frac{b^{2}}{4p^{2}}%
\right)
\end{equation}
with $\nu=0,p=1/2$ to express (by changing the order of integration)
\begin{align}
I(x)&=\int_0^\infty k e^{-k^2} J_0(kx)Y_0(kx)~\mathrm{d}k \\
&=\frac{1}{2}\int_0^\infty te^{-t^2/4}\,dt\int_0^\infty k J_0(kt) J_0(kx)Y_0(kx)~\mathrm{d}k
\end{align}
The integral of the triple product of Bessel function is (for a proof see Watson, "Theory of Bessel functions", §13.46):
\begin{equation}
\int_{0}^{\infty}Y_{\nu}\left(ak\right)J_{\nu}\left(bk\right)J_{\nu}\left(ck%
\right)k^{1+\nu}\mathrm{d}k=\begin{cases}-\dfrac{(abc)^{\nu}(-A)^{-\nu-\frac{1%
}{2}}}{\pi^{\frac{1}{2}}2^{\nu+1}\Gamma\left(\frac{1}{2}-\nu\right)},&0<a<|b-c%
|,\\
0,&|b-c|<a<b+c,\\
\dfrac{(abc)^{\nu}(-A)^{-\nu-\frac{1}{2}}}{\pi^{\frac{1}{2}}2^{\nu+1}\Gamma%
\left(\frac{1}{2}-\nu\right)},&a>b+c.\end{cases}
\end{equation} where $A=s(s-a)(s-b)(s-c), s=(a+b+c)/2$. Here, we choose $\nu=0,a=x,b=t,c=x$, and thus $A=-\left( t^2/4-x^2 \right)t^2/4$. This integral is non zero iff $t>2x$. In this case,
\begin{equation}
\int_{0}^{\infty}Y_{\nu}\left(ak\right)J_{\nu}\left(bk\right)J_{\nu}\left(ck%
\right)k^{1+\nu}\mathrm{d}k=-\frac{t^{-1}\left( t^2/4 -x^2\right)^{-1/2}}{\pi}
\end{equation}
Using simple changes of variables and an integral representation of the modified Bessel function (DLMF), we have then
\begin{align}
I(x)&=-\frac{1}{2\pi}\int_{2x}^\infty \frac{e^{-t^2/4}}{\sqrt{t^2/4-x^2}}\,dt\\
&=\frac{-1}{\pi}\int_0^\infty e^{-x^2\cosh^2u}du\\
&=-\frac{e^{-x^2/2}}{\pi}\int_0^\infty e^{-\frac{x^2}{2}\cosh 2u}\,du\\
&=-\frac{e^{-x^2/2}}{2\pi}\int_0^\infty e^{-\frac{x^2}{2}\cosh v}\,dv\\
&=-\frac{e^{-x^2/2}}{2\pi}K_0\left( \frac{x^2}{2} \right)
\end{align}
as expected.
