For simplicity assume $f$ is smooth. It turns out you can always approximate a subharmonic function by smooth subharmonic functions, so you aren't really losing any generality by doing so.
$(iv)\iff(iii)$: This is more or less by definition. Suppose $\phi\in C_c^\infty(U)$ is such that $\phi\geq 0$. Integration by parts says $$\int_Uf\Delta\phi = \int_U\phi\Delta f.$$ Suppose $(iv)$ holds, i.e., $\Delta f\geq 0$. Since $\phi\geq 0$, the right hand side of the equality is $\geq 0$, hence so is the left. Now assume $(iii)$. We know the left hand side is always $\geq 0$, and hence $\int \phi\Delta f\geq 0$ whenever $\phi\geq 0$ is compactly supported smooth. This means $\Delta f\geq 0$, as otherwise you could choose an appropriately constructed $\phi$ for which the right hand side would be negative.
You say you have shown $(i)\iff (iv)$, so this shows $(i)\iff (iv)\iff (iii)$.
The proof of $(ii)\implies (i)$ that I know uses a bit more -- the solution of the Dirichlet problem. Suppose $\mathbb{B}_r(x)\subset U$, and let $u$ be a harmonic function on $\mathbb{B}_r(x)$ that agrees with $f$ on the boundary of $\mathbb{B}_r(x)$. Because of $(ii)$, one has $f\leq u$ on all of $\mathbb{B}_r(x)$. In particular, $$f(x)\leq u(x) = \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x)} u(y)\,dS(y) = \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x)} f(y)\,dS(y).$$
It only remains to show that something implies $(ii)$. We'll do $(i)\implies (ii)$. Suppose that $(ii)$ is false. Let $D\Subset U$ be a domain and $u\colon D\to \mathbb{R}$ a harmonic function on $D$ which agrees with $f$ on $\partial D$, but for which there is a point $x\in D$ such that $f(x)>u(x)$. Let $V$ be the collection of points $\{x\in D : f(x)>u(x)$}. This is an open set which is nonempty by assumption. Fix an $x_0\in V$ for which $f(x_0) - u(x_0)$ is maximal. Then for small $r>0$, one has that $\mathbb{B}_r(x_0)\subset V$, and thus $$f(x_0) - u(x_0)\leq \frac{1}{n\alpha(n)r^{n-1}}\int_{\partial\mathbb{B}_r(x_0)} f(y) - u(y)\,dS(y).$$ On the other hand, the integrand $f(y)-u(y)$ is $\leq f(x_0) - u(x_0)$ by maximality, so the only way for this inequality to hold is if $f(y) - u(y)\equiv f(x_0) - u(x_0)$ for $y\in\partial\mathbb{B}_r(x_0)$. Since $r>0$ could have been anything (small), it follows that $f(x) - u(x)$ is constant in a neighborhood of $x_0$. By a connectedness argument, we see that $f(x) - u(x)$ is constant in $D$. However, $f$ and $u$ agree on $\partial D$, so it must be that $f\equiv u$ on $D$, contradicting our assumption that there is some $x$ with $f(x)>u(x)$. We conclude that $(i)\implies (ii)$.
No. A convex function is continuous in the interior of its domain, but it need not be continuous on the boundary.
For example, with $A = \{ x \in \mathbb{R}^n : \lVert x\rVert \leqslant 1\}$, where $\lVert \cdot\rVert$ is the Euclidean norm (or any strictly convex norm), the function
$$f(x) = \begin{cases}0 &, \lVert x\rVert < 1\\ g(x) &, \lVert x\rVert = 1 \end{cases}$$
is convex for every $g \colon \partial A \to [0,\infty)$.
Best Answer
Using that $u$ is midpoint-convex works in higher dimensions as well.
$y \mapsto x - (y-x) = 2x-y$ maps the sphere $\partial B_r(x)$ bijectively onto itself (each point is mapped to the “opposite” point on the sphere). It follows that $$ \int_{\partial B_r(x)} u(y) \, dy = \int_{\partial B_r(x)} u(2x-y) \, dy $$ and therefore $$ \int_{\partial B_r(x)} u(y) \, dy = \int_{\partial B_r(x)} \frac 12\bigl(u(y) + u(2x-y)\bigr) \, dy \\ \ge \int_{\partial B_r(x)} u\left(\frac{y + (2x-y)}{2}\right) \, dy = \int_{\partial B_r(x)} u(x) \, dy = |\partial B_r(x)| \cdot u(x) \, . $$