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.)