Define
$$
\newcommand{\dashint}{\mathchoice{\rlap{\,\,-}\int}{\rlap{\,-}\int}{\rlap{\ -}\int}{\rlap{\,-}\int}}
\dashint_Af(x)\,\mathrm{d}x=\frac1{|A|}\int_Af(x)\,\mathrm{d}x
$$
Using $(1)$ and $(2)$ from this answer and $|S(r,p)|=\omega_{n-1}r^{n-1}$, we get that
$$
\begin{align}
&r_2^{n-1}\frac\partial{\partial r}\dashint_{S(r_2,p_2)}u(x)\,\mathrm{d}\sigma-r_1^{n-1}\frac\partial{\partial r}\dashint_{S(r_1,p_1)}u(x)\,\mathrm{d}\sigma\\
&=\frac1{\omega_{n-1}}\int_{B(r_2,p_2)\setminus B(r_1,p_1)}\Delta u(x)\,\mathrm{d}x\\[6pt]
&=0\tag1
\end{align}
$$
where $r_k\gt|p_k|$ (the origin is inside both spheres).
Thus, for some constant $C$, independent of $p$ (as long as $r\gt|p|$),
$$
r^{n-1}\frac\partial{\partial r}\dashint_{S(r,p)}u(x)\,\mathrm{d}\sigma=C\tag2
$$
Therefore,
$$
\dashint_{S(r,p)}u(x)\,\mathrm{d}\sigma
=\left\{\begin{array}{}
A(p)-\frac{C}{(n-2)\,r^{n-2}}&\text{if }n\ge3\\
A(p)+C\log(r)&\text{if }n=2
\end{array}\right.
\tag3
$$
If $u$ is bounded, then by considering $p=0$, we get $C=0$, and therefore, for $r\gt|p|$,
$$
\dashint_{S(r,p)}u(x)\,\mathrm{d}\sigma=A(p)\tag4
$$
Since $u$ is harmonic away from $0$, the Mean Value Property says that for $r\lt|p|$,
$$
\dashint_{S(r,p)}u(x)\,\mathrm{d}\sigma=u(p)\tag5
$$
However, because $u$ is smooth away from $0$, and bounded, $\dashint_{S(r,p)}u(x)\,\mathrm{d}\sigma$ is a continuous function of $r$. That is, $A(p)=u(p)$ and we have $(5)$ for all $r$.