$\Delta \mathbf n = -2 \mathbf n$ on the Euclidean sphere Let us consider the Euclidean two-sphere, defined by the embedding in the three dimensional Euclidean space as
$$
\mathbf n \cdot \mathbf n = 1\,,
$$
where $\cdot$ denotes the standard scalar product. The metric on the sphere, in some coordinates $x^i$, is expressed as 
$$
\gamma_{ij}=\mathbf e_i \cdot \mathbf e_j\,,
$$
where $\mathbf e_i=\partial_i\mathbf n$. For instance, in the standard spherical coordinates 
$$
\mathbf n=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)
$$
and
$$
\gamma_{\theta\theta}=1\,,\qquad \gamma_{\phi\phi}=\sin^2\theta\,,\qquad
\gamma_{\theta\phi}=0\,.
$$
We define the Laplace-Beltrami operator on the sphere by
$
\Delta = \gamma^{ij} D_iD_i
$,
where $\gamma^{ij}$ is its inverse and $D_i$ is the associated Levi-Civita connection.
I would like to prove that
$$\Delta \mathbf n = -2 \mathbf n$$
and that in higher dimensions the same holds with $2$ replaced by the dimension of the sphere. I came to believe that this is true by an explicit check in spherical coordinates in dimensions $3$, $4$ and $6$.
Considering that parallel transport of a given tangent vector $\mathbf v$ defined at the point $x+dx$ to the point $x$ is defined by keeping constant its embedding components and then projecting it on the sphere at the point $x$, we have
$$
\mathbf v_{\parallel}(x+dx,x)=\mathbf v(x+dx)-\mathbf v(x+dx)\cdot \mathbf n (x)\, \mathbf n(x)
$$
hence
$$
D_i\mathbf v\, dx^i = \mathbf v_{\parallel}(x+dx,x)- \mathbf v (x)= (\partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v\, \mathbf n)dx^i
$$
where we have used $\mathbf n \cdot \partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v=0$ and
$$
D_i\mathbf v = \partial_i\mathbf v+\partial_i\mathbf n \cdot \mathbf v\, \mathbf n\,.
$$
Applying this to the basis vectors $\mathbf e_j =\partial_j\mathbf n$ affords
$$
D_i\mathbf e_j = \partial_i \mathbf e_j+\gamma_{ij}\mathbf n\,.
$$
But unfortunataly I am not able to go further.

You can also develop in spherical harmonics. This works equally well on spheres of arbitrary dimension. In Cartesian coordinates $x_1, x_2, \ldots, x_d$, the vector $\boldsymbol{n}$ on $\mathbb{S}^{d-1}$ is
$$
\boldsymbol n = (x_1, x_2, \ldots, x_d), $$
and since each entry is, manifestly, a homogeneous harmonic polynomial of degree $1$, $\boldsymbol n$ is a spherical harmonic of degree $1$. Now, it is well-known that a spherical harmonic of degree $\ell$ is an eigenvector of $\Delta_{\mathbb{S}^{d-1}}$ with eigenvalue $-\ell(\ell+d-2)$. Specializing to $\ell=1$ completes the proof that
$$
\Delta_{\mathbb{S}^{d-1}} \boldsymbol n= -(d-1)\boldsymbol n.$$
REMARK. I will prove that the eigenfunctions of the Laplace-Beltrami are precisely the homogeneous harmonic polynomials, as stated above. This illustrates a computational technique that can be often used in practice for the Laplace-Beltrami operator.
Let $H_\ell=H_\ell(x_1, x_2, \ldots, x_d)$ be a homogeneous harmonic polynomial of degree $\ell$. Then, by harmonicity, $\Delta_{\mathbb R^d}H_\ell = 0$. Expanding the Laplacian in polar coordinates yields
$$\tag{*}\Delta_{\mathbb R^n} H_{\ell}=
\frac{1}{r^{d-1}}\partial_r(r^{d-1}\partial_r H_\ell) + \frac{1}{r^2}\Delta_{\mathbb S^{d-1}}H_\ell=0.$$
Since $H_\ell$ is homogeneous, by the Euler theorem on homogeneous functions we have that $x\cdot \nabla H_\ell=\ell H_\ell$, that is
$$
r\partial_r H_\ell = \ell H_\ell.$$
Plugging this into (*), we arrive at
$$
\frac{\ell(\ell+d-2)}{r^2}H_\ell + \frac{\Delta_{\mathbb{S}^{d-1}}H_\ell}{r^2}=0, $$
from which our claim immediately follows.
