Yes.

It follows from the spectral theorem that for any $h \in L^2$, we have $P_t h \in D(\Delta)$.

I think this is easiest to see from the multiplication operator form of the spectral theorem; i.e., there is a measure space $(\hat{X}, \mu)$, a measurable function $g : \hat{X} \to \mathbb{R}$, and a unitary operator $U : L^2(X) \to L^2(\hat{X})$ such that $\Delta = U^{-1} M_g U$. Here $M_g$ is the multiplication operator on $L^2(\hat{X})$ defined by $M_g v = gv$ with domain $D(M_g) = \{ v \in L^2(\hat{X}) : gv \in L^2(\hat{X})\}$. Since $\Delta$ is a negative definite operator, $g \le 0$ almost everywhere.

Now it is easy to check that $P_t = U^{-1} M_{e^{tg}} U$. But notice that $g e^{tg}$ is a bounded function. Thus for any $v \in L^2(\hat{X})$, we have $M_{e^{tg}} v \in D(M_g)$. Using $U$ to move back to $L^2(X)$, we see that for any $h \in L^2(X)$, $P_t h \in D(\Delta)$.

A similar argument shows that $P_t h \in D(\Delta^\infty)$.

To show $P_t h \in L^\infty$, one could use the heat kernel bounds shown by Li and Yau:

Li, Peter; Yau, Shing-Tung.
On the parabolic kernel of the Schrödinger operator.
Acta Math. 156 (1986), no. 3-4, 153–201.

They show, among other things, that in a complete Riemannian manifold with Ricci curvature bounded from below, we have Gaussian upper estimates for the heat kernel. In particular, for each $t$, $p_t$ is bounded. From this it is obvious that $P_t$ maps $L^1$ into $L^\infty$.

Now $P_t$ also maps $L^\infty$ into $L^\infty$, so it maps $L^1 + L^\infty$ into $L^\infty$; in particular it maps $L^2$ into $L^\infty$.

Since $P_t$ maps $L^\infty$ into $L^\infty$ and is symmetric, it also maps $L^1$ into $L^1$. So it maps $L^1$ into $L^1 \cap L^\infty$; in particular it maps $L^1$ into $L^2$. In particular, since $p_t(x_0, \cdot) = P_{t/2} [p_{t/2}(x_0, \cdot)]$, where $p_{t/2}(x_0, \cdot) \in L^1$, we have $p_t(x_0, \cdot) \in L^2$. Then by using $p_t(x_0, \cdot) = P_{t/2} [p_{t/2}(x_0, \cdot)]$ again, we get $p_t(x_0, \cdot) \in D(\Delta)$.

Finally, we note that for $h \in D(\Delta)$, we have $\Delta P_t h = P_t \Delta h$. So we have
$$\Delta [p_t(x_0, \cdot)] = \Delta P_{t/2} [p_{t/2}(x_0, \cdot)] = P_{t/2} \Delta [p_{t/2}(x_0, \cdot)].$$
Since $p_{t/2}(x_0, \cdot) \in D(\Delta)$ as previously argued, we have $\Delta p_{t/2}(x_0, \cdot) \in L^2$. Since $P_{t/2}$ maps $L^2$ into $L^\infty \cap D(\Delta)$, we have $\Delta [p_t(x_0, \cdot)] \in L^\infty \cap D(\Delta)$ as desired.