Positive harmonic function on $\mathbb{R}^n$ is a constant? Is it true that a positive harmonic function on $\mathbb{R}^n$ must be a constant? How might we show this? The mean value property seems not to be the way...for that we would need boundedness.

Since it is not clear whether the Wikipedia proof uses boundedness or not, please allow me to give a detailed proof that only uses nonnegativity.
Let $u$ be a nonnegative harmonic function in $\mathbb{R}^n$, and let $x,y\in\mathbb{R}^n$. Denote by $B_R(y)$ the ball of radius $R>0$ centred at $y$, and similarly by $B_r(x)$ the ball of radius $r>0$ centred at $x$. We consider such balls with $R>0$ large and $r=R-|x-y|$, that is, $B_r(x)\subset B_R(y)$ and the boundary of $B_r(x)$ is tangent to the boundary of $B_R(y)$.
Then we have
$$
u(x) = \frac1{|B_r|}\int_{B_r(x)}u \leq \frac1{|B_r|}\int_{B_R(y)}u = \frac{|B_R|}{|B_r|}u(y),
$$
where we have used the mean value property in the first and last equalities, and the nonnegativity of $u$ in the second inequality. Now we send $R\to\infty$ and get 
$$
u(x) \leq u(y).
$$
As $x$ and $y$ were arbitrary, we conclude that $u$ is constant.
Note that the Harnack inequality is lurking behind this proof. Also, it was crucial to consider balls with different radii in order to make use of nonnegativity. The Wikipedia proof uses balls of equal radii.
