It was requested by another user that I post the method to show that $$\lim_{t \to 0}u(t,x)=u_0(x)$$ assuming that $u_0$ has a derivative, not necessarily continuous, which is bounded on every finite interval. I'll also asume that
$$
\int_{-\infty}^{\infty} |u_0(x)|\,dx < \infty.
$$
My apologies if you don't find this relevant.
Fix $x$. Let $\lambda = 1/(4t)$ and make the change of variables $s = y-x$. For convenience, define $v_x(s) = u_0(x+s)$. This gives
$$
\begin{align*}
\int_{-\infty}^{\infty} e^{-(x-y)^2/(4t)} u_0(y)\,dy &= \int_{-\infty}^{\infty} e^{-\lambda s^2} v_x(s) \,ds \\
&= \int_{-\infty}^{\infty} e^{-\lambda s^2} \Bigl[v_x^{\text{odd}}(s) + v_x^{\text{even}}(s)\Bigr] \,ds \\
&= 2\int_{0}^{\infty} e^{-\lambda s^2} v_x^{\text{even}}(s) \,ds,
\tag{1}
\end{align*}
$$
where
$$
v_x^{\text{even}}(s) \triangleq \frac{v_x(s)+v_x(-s)}{2} \qquad \text{and} \qquad v_x^{\text{odd}}(s) \triangleq \frac{v_x(s)-v_x(-s)}{2}.
$$
Making the change of variables $s = \sqrt{r}$ in $(1)$ gives
$$
\begin{align*}
2\int_{0}^{\infty} e^{-\lambda s^2} v_x^{\text{even}}(s) \,ds &= \int_{0}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr \\
&= \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr + \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr,
\qquad \tag{2}
\end{align*}
$$
where $0 < \delta < \infty$ is fixed. We see that the right integral in $(2)$ is exponentially small as $\lambda \to \infty$ by the calculation
$$
\left|\int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr\right| \leq e^{-\delta \lambda} \int_{\delta}^{\infty} r^{-1/2} \left|v_x^{\text{even}}\left(\sqrt{r}\right)\right| \,dr.
\tag{3}
$$
The finiteness of this last integral follows from the assumption that $u_0$ is integrable.
Now, by assumption we can appeal to Taylor's theorem with remainder to write
$$
\begin{align*}
v_x(s) &= v_x(0) + R_x(s) \\
&= u_0(x) + R_x(s),
\end{align*}
$$
where
$$
|R_x(s)| \leq \sup_{|\tau| < \delta} |u_0'(x+\tau)| \cdot |s|
$$
for $|s| \leq \delta$, so that
$$
v_x^{\text{even}}(s) = u_0(x) + \hat{R}_x(s),
\tag{4}
$$
where
$$
|\hat{R}_x(s)| \leq \sup_{|\tau| < \delta} |u_0'(x+\tau)| \cdot |s|
$$
for $|s| \leq \delta$. Substituting $(3)$ into the left integral in $(2)$ we get
$$
\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} v_x^{\text{even}}\left(\sqrt{r}\right) \,dr = u_0(x) \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \,dr + \int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \hat{R}_x\left(\sqrt{r}\right) \,dr.
\quad \tag{5}
$$
We can bound the integral on the right,
$$
\left|\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \hat{R}_x\left(\sqrt{r}\right) \,dr\right| \leq \sup_{|\tau| < \delta} |u_0'(x-\tau)| \int_0^\delta e^{-\lambda r}\,dr = O(1/\lambda),
\tag{6}
$$
and the integral on the left can be rewritten as
$$
\begin{align*}
\int_{0}^{\delta} e^{-\lambda r} r^{-1/2} \,dr &= \int_{0}^{\infty} e^{-\lambda r} r^{-1/2} \,dr - \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr \\
&= \sqrt{\frac{\pi}{\lambda}} - \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr,
\tag{7}
\end{align*}
$$
where
$$
0 \leq \int_{\delta}^{\infty} e^{-\lambda r} r^{-1/2} \,dr \leq \frac{1}{\sqrt{\delta}\lambda} e^{-\delta\lambda}
\tag{8}
$$
is exponentially small as $\lambda \to \infty$.
Finally, combining the estimates $(3)$, $(6)$, and $(8)$ with the expressions $(1)$, $(2)$, $(5)$, and $(7)$, we get
$$
\int_{-\infty}^{\infty} e^{-\lambda\,(x-y)^2} u_0(y)\,dy = \sqrt{\frac{\pi}{\lambda}}\,u_0(x) + O\left(\frac{1}{\lambda}\right)
$$
as $\lambda \to \infty$. Equivalently,
$$
\frac{1}{\sqrt{4\pi t}} \int_{-\infty}^{\infty} e^{-(x-y)^2/(4t)} u_0(y)\,dy = u_0(x) + O\left(\sqrt{t}\right)
$$
as $t \to 0^+$.
(The result is surely true for lighter assumptions on $u_0$, perhaps if $|u_0(x)| \leq a e^{b x^2}$. This could add to the calculations considerably and I wanted to keep things more simple than complicated here.)