1) By first removing the cut locus of $p$, the remaining subset of $M$ is an open subset diffeomorphic to some open subset of $T_p M$. Furthermore, if you place yourself in spherical normal coordinates $(r, \sigma)$ (where $\sigma$ encodes all the angular variables), the function $d_p ^2$ becomes $r^2$, which is obviously smoooth.
For a rigorous statement (but following the same idea), see theorem 3.6 at page 166, in chapter IV of "Foundations of Differential Geometry", volume 1, by Kobayashi and Nomizu.
2) In the same spherical normal coordinates, the Laplacian looks like
$$\Delta f \ (q) = \frac {\partial ^2 f} {\partial r^2} (q) + H(p,q) \frac {\partial f} {\partial r} (q) + \frac 1 {r^2} \Delta_S f \ (q) ,$$
where $r = d(p,q)$, $S$ is the geodesic sphere of radius $r$ centered at $p$, $\Delta_S$ is the Laplacian restricted to $S$ and $H(p,q)$ is the mean curvature at $q$ of $S$. It follows that
$$\Delta d_p ^2 \ (q) = \Delta r^2 = 2 + H(p,q) 2r = 2 + 2 H(p,q) d_p (q) .$$
In $\Bbb R^n$, the mean curvature of the sphere of radius $r$ is $H = \frac {n-1} r$, so the above becomes $2n$. Otherwise, as your intuition hinted at, the curvature does play a role and makes the result non-constant, in general.
For an alternative formulation of the Laplacian see also page 48 of "Spectral Theory and Geometry" edited by E.B. Davies and Y. Safarov.
If you use normal coordinate, the computation along this line will be easier: here we have $d=\sqrt{x_1^2+...+x_n^2}$, $\Gamma_{ij}^k(p)=0$ so $\Gamma_{ij}^k=O(d)$; and $\frac{\partial g_{ij}}{\partial x_k}(p)=0$, so $g_{ij}=\delta_{ij}+O(d^2)$. Take inverse we get also $g^{ij}=\delta_{ij}+O(d^2)$. Now from
$\Delta d=g^{ij}(d_{ij}-\Gamma_{ij}^kd_k)$, by direct calculation we note $d_k=O(1)$, so $\Gamma_{ij}^kd_k\to 0$. Also by direct computation we see $d_{ij}=O(d^{-1})$; so $g^{ij}d_{ij}=\delta_{ij}d_{ij}+O(d^2)O(d^{-1})$, with $ O(d^2)O(d^{-1})\to 0$. Thus at distance $d$, $\Delta d$ differs from the Euclidean version by $O(d)\to 0$, the result follows from the Euclidean case.
Best Answer
Isn't it immediate from the formula for Laplacian in a normal ball?
We have $$ \Delta_g = -\partial_r^2+\left(\frac{n-1}{r}+\frac{\partial_r J}{J}\right)\partial_r+\Delta_{S_g^{n-1}(r)} $$ where $J=\det(d\exp_p)$. Hitting this on $r$, we get $$ \Delta_g r=-\frac{n-1}{r}-\frac{\partial_r J}{J} $$ But $J\to 1$ as $r\downarrow 0$ so there is no chance of cancelling the $(n-1)/r$ terms on the right.