In fact, Poisson integral formula implies the following lemma.
(Harnack’s Inequality) Suppose $u$ is harmonic in $B_R(x_0)$ and $u\geq0$.Then there holds
$$\left(\frac{R}{R+r}\right)^{n-2}\frac{R-r}{R+r}u(x_0)\leq u(x)\leq \left(\frac{R}{R-r}\right)^{n-2}\frac{R+r}{R-r}u(x_0) $$
where $r=|x-x_0|<R.$
Proof. Without loss of a generality, we may assume $x_0=0$ and $u\in C(\overline{B}_R).$ By the Poisson integral formula, we have
$$u(x)=\frac{R^2-|x|^2}{n\alpha(n)R}\int_{\partial B_R}\frac{u(y)}{|x-y|^n}\mathrm{~d}S$$
where $\alpha(n)$ is the volume of unit sphere.
Since $R-r\leq |x-y|\leq R+r$ with $|y|=R$, we have
\begin{align*}
\frac{1}{n\alpha(n)R}\cdot\frac{R-|x|}{R+|x|}&\left(\frac{1}{R+|x|}\right)^{n-2}\int_{\partial B_R} u(y)\mathrm{~d}S\leq u(x)\\&\leq \frac{1}{n\alpha(n)R}\cdot\frac{R+|x|}{R-|x|}\left(\frac{1}{R-|x|}\right)^{n-2} \int_{\partial B_R} u(y)\mathrm{~d}S.
\end{align*}
Then the
Mean Value Property gives us that
$$u(0)=\frac{1}{n\alpha(n)R^{n-1}}\int_{\partial B_R}u(y)\mathrm{~d}S,$$
which completes the proof.
Hence Liouville theorem is a corollary of the lemma.
In particular, We assume $u\geq 0$ in $\mathbb{R^n}$. Taking any point $x\in\mathbb{R^n}$ and applying the lemma to any ball $B_R(0)$ with $R > |x|$, we obtain
$$\left(\frac{R}{R+r}\right)^{n-2}\frac{R-r}{R+r}u(0)\leq u(x)\leq \left(\frac{R}{R-r}\right)^{n-2}\frac{R+r}{R-r}u(0), $$
which yields $u(x)=u(0)$ by letting $R\to +\infty$.
If we don’t use Poisson's formula, we shall give another proof. The Mean Value Formula implies the following lemma.
Suppose $u\in C(\overline{B}_R(x_0))$ is a nonnegative harmonic function in $B_R=B_R(x_0)$. Then there holds
$$|Du(x_0)|\leq \frac{n}{R}u(x_0).$$
Proof. Since u is smooth in $B_R$, we know $\Delta (D_{x_i}u)=0$, that is $D_{x_i}u$ is also harmonic in $B_R$. Hence $D_{x_i}u$ satisfies mean value formula. Then by divergence theorem and the nonnegativeness of $u$ we have
\begin{align}|D_{x_i}u(x_0)|=\left|\frac{1}{\alpha(n)R^n}\int_{B_R}D_{x_i}u(y)\,dy\right|&=\left|\frac{1}{\alpha(n)R^n}\int_{\partial B_R} u(y) v_{x_i} \,dS\right|\\&\leq \frac{1}{\alpha(n)R^n}\int_{\partial B_R} u(y) \,dS\\&=\frac{n}{R}u(x_0)
\end{align}
where in the last equality we used the mean value property.
Then letting $R\to +\infty$, we get $Du(x)=0$ for all $x\in\mathbb{R^n},$ which implies the Liouville theorem.
All above are from Qin Han and Fanghua Lin’s elliptic partial differential equations.