Harmonic in punctured ball and bounded implies harmonic in ball.

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

The function $u$ is harmonic on a neighborhood of $0$ with $u(x) |x|^{n-2} \to 0$ as $x\to 0$, so it extends to a harmonic function across $0$. (This is a standard result about harmonic functions, but the proof is basically solving the Dirichlet problem on a small punctured ball around $0$ with boundary value $u$, then using the maximum principle and the condition above to ensure that this new solution must just be $u$ itself.)