Yes, the rest mass of a system such as a hydrogen atom can be calculated in a fully gauge-invariant way, by integrating the $T^{00}$ component of its gauge-invariant Hilbert energy-momentum-stress tensor $T^{\mu\nu}$. This component is the energy density, and integrating it over all space gives the total gauge-invariant energy of the system.
The "canonical" energy-momentum-stress tensor' produced by Noether's Theorem is guaranteed to be conserved if the action for the system is invariant under temporal and spatial translations, but it isn't necessarily gauge-invariant or even symmetric. Thus it, and the energy found by integrating it, has no physical significance.
However, there is a well-known procedure for producing physically meaningful energy-momentum-stress tensors which are not only conserved but also manifestly gauge-invariant, manifestly Lorentz-covariant, and manifestly symmetric. These tensors do have physical significance: they represent the local energy density, momentum density, etc. of matter.
From the appearance of this kind of $T^{\mu\nu}$ on the right-hand-side of Einstein's field equations for General Relativity, it is obvious that energy density, momentum density, etc. must be measurable, because they create curvature in spacetime. Curvature is measurable, and therefore $T^{\mu\nu}$ is measurable. Thus is it not the case that only differences in energies are physically meaningful. The "absolute" value of the energy density is meaningful, because it curves spacetime.
There is a standard procedure due to Hilbert and Einstein for finding stress-tensors that are conserved, gauge-invariant, Lorentz-covariant, and symmetric, namely by functional integration of the Lagrangian density with respect to the metric:
$$T_{\mu\nu} = \frac{-2}{\sqrt{-g}} \frac{{\delta\mathcal{L}_{\text{matter}}\sqrt{-g}}}{\delta g^{\mu\nu}}.$$
When one uses such a $T^{\mu\nu}$ for a system of point particles and electromagnetic fields, the result for the hydrogen atom (in the order-$\alpha^2$ approximation that the question uses) is identical to the "naïve" calculation presented in the question. But it is completely gauge-invariant from beginning to end.
What one finds is that $U(r) = q_1 q_2/r$, previously considered the "electrostatic potential energy of the proton and electron" is in fact the gauge-invariant position-dependent part of the rest energy residing in the electrostatic field of the proton and electron. The fact that it looks like the product of one charge with the gauge-dependent electrostatic potential of the other charge is irrelevant. The calculation below does not involve any potentials at all.
There is also an infinite, position-independent, gauge-invariant, constant part of the electrostatic field energy, which simply renormalizes the masses of the two particles.
Thus electrostatic "potential energy" can be understood as the position-dependent part of the gauge-invariant rest energy residing in the electrostatic field. It really goes to zero at infinity, because a gauge-invariant calculation tells us so. Going to zero at infinity is not simply a convention.
If you want to say that the mass of a proton and an electron are their standard measured values, then you must take the "potential energy" between them to be exactly $-e^2/r$ and not add a constant or any other gauge-dependent term.
Here are the mathematical details:
Consider a collection of point particles with masses $m_i$ and charges $q_i$, moving under the influence of an electromagnetic field, which could be a combination of the fields of the particles themselves and some external field.
The conserved, symmetric, manifestly-covariant, and gauge-invariant energy-momentum-stress tensor for a system of point particles and electromagnetic fields divides cleanly into two terms,
$$T^{\mu\nu} = T_{(p)}^{\mu\nu} + T_{(f)}^{\mu\nu},$$
where the particles-only term is
$$T_{(p)}^{\mu\nu} = \sum_i m_i \int \dot{X}_i^\mu (\tau) \dot{X}_i^\nu (\tau) \; \delta^{(4)}(x-X_i(\tau)) \; d\tau$$
and the fields-only term is
$$T_{(f)}^{\mu\nu}=\frac{1}{4\pi} \left( F^{\mu\alpha} {F^\nu}_\alpha - \frac{1}{4} \eta^{\mu\nu} F^{\alpha\beta} F_{\alpha\beta} \right).$$
Here $m_i$ is the mass of the $i^\text{th}$ particle, $X_i(\tau)$ is its world-line as a function of the proper time $\tau$ along the worldline, and $F^{\mu\nu}$ is the gauge-invariant electromagnetic field tensor. Note that the gauge-dependent electromagnetic potential appears nowhere in this energy-momentum-stress tensor.
The $T^{00}$ component for the particles can be written, after doing the integration over $\tau$, as
$$T_{(p)}^{00} = \sum_i m_i( 1-\mathbf{v}_i^2)^{-1/2} \; \delta^{(3)}(\mathbf{x}-\mathbf{X}_i(t)),$$
and integrating this over space gives the energy of the particles,
$$E_{(p)} = \sum_i m_i( 1-\mathbf{v}_i^2)^{-1/2} = \sum_i m_i + \sum_i \frac{1}{2}m_i \mathbf{v}^2 + \dots \; .$$
Obviously this is the usual expansion into the energy of the rest masses plus the kinetic energy. We're doing a non-relativistic approximation and don't need the higher terms. In the case of the hydrogen atom, the energy of the particles in the center-of-mass frame takes the form
$$E_{(p)} = m_p + m_e + \frac{\mathbf{p}^2}{2\mu}$$
where $\mu = m_p m_e/(m_p + m_e)$ is the reduced mass and $\mathbf{p}$ is the relative momentum. The expectation value of this third term is what was denoted $\langle K \rangle$ in the question. So we have reproduced the first three terms for the rest energy of hydrogen.
The $T^{00}$ component for the fields can be written in terms of the electric field $\mathbf{E}$ and the magnetic field $\mathbf{B}$,
$$T_{(f)}^{00} = \frac{\mathbf{E}^2 + \mathbf{B}^2}{8 \pi},$$
and integrating this over space gives the energy of the fields,
$$E_{(f)} = \int \frac{\mathbf{E}^2 + \mathbf{B}^2}{8 \pi} \, dV.$$
For the order-$\alpha^2$ approximation of the hydrogen rest energy that we care about, we can calculate this field energy by ignoring the motion of the charges. The electric field is just the usual Coulomb one for a static charge, and there is no magnetic field. There is no external field to consider.
The fields of $q_1$ are
$$\mathbf{E}_1 = q_1 \frac{\mathbf{r}-\mathbf{r_1}}{|\mathbf{r}-\mathbf{r_1}|^3}, \quad \mathbf{B}_1 = 0,$$
and the fields of $q_2$ are
$$\mathbf{E}_2 = q_2 \frac{\mathbf{r}-\mathbf{r_2}}{|\mathbf{r}-\mathbf{r_2}|^3}, \quad \mathbf{B}_2 = 0.$$
The field energy obviously breaks into three integrals:
$$E_f = \int \frac{(\mathbf{E}_1 + \mathbf{E}_2)^2}{8 \pi} \, dV
= \int \frac{\mathbf{E}_1^2}{8 \pi} \, dV + \int \frac{\mathbf{E}_1 \cdot \mathbf{E}_2}{4 \pi} \, dV + \int \frac{\mathbf{E}_2^2}{8 \pi} \, dV.$$
The first integral,
$$E_{f,1} = \int \frac{\mathbf{E}_1^2}{8 \pi} \, dV = \frac{q_1^2}{8\pi} \int \frac{1}{|\mathbf{r}-\mathbf{r_1}|^4} \, dV = \frac{q_1^2}{2} \int_0^\infty \frac{dr}{r^2},$$
diverges. It is the classical electrostatic self-energy of point charge $q_1$ interacting with its own field. It has nothing to do with $q_2$ and does not depend on the distance between the charges or the position of $q_1$. It is rest energy that is just as intrinsic to $q_1$ as its mass-energy is, and this energy simply renormalizes the mass $m_1$, in the same way as happens in QED.
The third integral,
$$E_{f,2} = \int \frac{\mathbf{E}_2^2}{8 \pi} \, dV = \frac{q_2^2}{8\pi} \int \frac{1}{|\mathbf{r}-\mathbf{r_2}|^4} \, dV = \frac{q_2^2}{2} \int_0^\infty \frac{dr}{r^2},$$
similarly diverges and simply renormalizes $m_2$.
The second integral,
$$E_{f,12} = \int \frac{\mathbf{E}_1 \cdot \mathbf{E}_2}{4 \pi} \, dV = \frac{q_1 q_2}{4\pi} \int \frac{\mathbf{r}-\mathbf{r_1}}{|\mathbf{r}-\mathbf{r_1}|^3} \cdot \frac{\mathbf{r}-\mathbf{r_2}}{|\mathbf{r}-\mathbf{r_2}|^3} \; d^3\mathbf{r} = \;?$$
is the interesting one. It involves both gauge-invariant fields $\mathbf{E}_1$ and $\mathbf{E}_2$ and can be described as a gauge-invariant position-dependent interaction energy. It looks like it might be divergent, but it turns out that this integral can be performed (see below) and the result is finite! In fact, it is just the usual "potential energy" of two point charges:
$$E_{f,12} = \frac{q_1 q_2}{|\mathbf{r_1}-\mathbf{r_2}|}$$
For the hydrogen atom, this means that the "potential energy" term $\langle -e^2/r \rangle$ is completely legitimate; it is the $r$-dependent part of the gauge-invariant field energy. The $r$-independent part of the gauge-invariant field energy is divergent and renormalizes the masses.
This term coming from the fields produces the fourth term in the hydrogen rest energy, which was previously denoted $\langle U\rangle$.
To do the integral, introduce Cartesian coordinates with $q_1$ is at $\mathbf{r}_1=a\hat{\mathbf{z}}$ and $q_2$ is at $\mathbf{r}_2=-a\hat{\mathbf{z}}$.
The electric field of $q_1$ is
$$\mathbf{E}_1 = q_1 \frac{x\hat{\mathbf{x}} + y\hat{\mathbf{y}} + (z-a)\hat{\mathbf{z}}}{[x^2 + y^2 + (z-a)^2]^{3/2}}$$
and the electric field of $q_2$ is
$$\mathbf{E}_2 = q_2 \frac{x\hat{\mathbf{x}} + y\hat{\mathbf{y}} + (z+a)\hat{\mathbf{z}}}{[x^2 + y^2 + (z+a)^2]^{3/2}},$$
so the field interaction energy is
$$E_{f,12} = \frac{q_1 q_2}{4 \pi} \int_{-\infty}^\infty dx \int_{-\infty}^\infty dy \int_{-\infty}^\infty dz \; \frac{x^2 + y^2 + z^2 - a^2}{[(x^2 + y^2 + z^2 + a^2)^2 - 4 a^2 z^2]^{3/2}}.$$
Converting to spherical polar coordinates, the integral becomes
$$E_{f,12} = \frac{q_1 q_2}{4 \pi} \int_0^\infty r^2 dr \int_0^\pi \sin\theta d\theta \int_0^{2\pi} d\phi \; \frac{r^2 - a^2}{[(r^2 + a^2)^2 - 4 a^2 r^2 \cos^2\theta]^{3/2}}.$$
The integral over $\phi$ just gives $2\pi$, and the integral over $\theta$ becomes elementary with the subsitution $u=\cos\theta$:
\begin{align*}
E_{f,12} &= \frac{q_1 q_2}{2} \int_0^\infty r^2 (r^2 - a^2) \; dr \int_{-1}^1 \frac{du}{[(r^2 + a^2)^2 - 4 a^2 r^2 u^2]^{3/2}} \\
&= \frac{q_1 q_2}{2} \int_0^\infty r^2 (r^2 - a^2) \; dr \left[ \frac{u}{(r^2 + a^2)^2 [(r^2 + a^2)^2 - 4 a^2 r^2 u^2]^{1/2}} \right]_{-1}^1 \\
&= q_1 q_2 \int_0^\infty \frac{r^2 dr}{(r^2 + a^2)^2} \frac{r^2 - a^2}{|r^2 - a^2|}.
\end{align*}
The integral over $r$ must be split into two parts, one from 0 to $a$, and one from $a$ to $\infty$.
$$E_{f,12} = q_1 q_2 \left( \int_0^a \frac{r^2 dr}{(r^2 + a^2)^2}(-1) + \int_a^\infty \frac{r^2 dr}{(r^2 + a^2)^2} (+1) \right).$$
The resulting integrals can be done with the substitution $r=a\tan{u}$ and give
\begin{align*}
E_{f,12} &= q_1 q_2 \left( \left[ \frac{r}{2(r^2 + a^2)} - \frac{\tan^{-1}(r/a)}{2 a} \right]_0^a + \left[ \frac{\tan^{-1}(r/a)}{2 a} - \frac{r}{2(r^2 + a^2)} \right ]_a^\infty \right) \\
&= q_1 q_2 \left[ \left( \frac{1}{4a} - \frac{\pi}{8a} \right) - (0) + \left( \frac{\pi}{4a} \right) - \left( \frac{\pi}{8a} - \frac{1}{4a} \right) \right] \\
&= \frac{q_1 q_2}{2a}
\end{align*}
Thus the final result is
$$E_{f,12} = \frac{q_1 q_2}{d}$$
where $d=2a$ is the separation between the charges.
The claimed result,
$$E_{f,12} = \int \frac{\mathbf{E}_1 \cdot \mathbf{E}_2}{4 \pi} \, dV = \frac{q_1 q_2}{4\pi} \int \frac{\mathbf{r}-\mathbf{r_1}}{|\mathbf{r}-\mathbf{r_1}|^3} \cdot \frac{\mathbf{r}-\mathbf{r_2}}{|\mathbf{r}-\mathbf{r_2}|^3} \; d^3\mathbf{r} = \frac{q_1 q_2}{|\mathbf{r_1}-\mathbf{r_2}|}$$
follows from the fact that this is a rotationally-invariant equation which we've verified with one particular choice of coordinates.
Addendum:
A much simpler approach is simply to realize that the Schrodinger equation is gauge-invariant. For an electron in an electromagnetic field described by a scalar potential $\phi(\mathbf{r})$ and a vector potential $\mathbf{A}(\mathbf{r})$, the Schrodinger equation is
$$\left[ \frac{1}{2m} \left( \frac{\hbar}{i}\nabla + \frac{e}{c} \mathbf{A}(\mathbf{r}) \right)^2 - e \phi(\mathbf{r}) \right] \psi(\mathbf{r}, t)=E \psi(\mathbf{r}, t).$$
When one makes the gauge transformations
$$\psi(\mathbf{r}, t) \rightarrow \exp({i \lambda(\mathbf{r}, t)}) \psi(\mathbf{r}, t)$$
$$\phi(\mathbf{r}) \rightarrow \phi(\mathbf{r}) + \frac{\hbar}{e} \frac{\partial\lambda(\mathbf{r}, t)}{\partial t}$$
$$\mathbf{A}(\mathbf{r}) \rightarrow \mathbf{A}(\mathbf{r}) - \frac{\hbar c}{e} \nabla \lambda(\mathbf{r}, t)$$
and uses the relationship
$$i \hbar \frac{\partial \psi(\mathbf{r}, t)}{\partial t} = E \psi(\mathbf{r}, t),$$
one finds that the equation is unchanged, showing that the energy $E$ is gauge-invariant. This is a standard part of most undergraduate courses in Quantum Mechanics.
Thus the original calculation in the question is completely gauge-invariant (although not manifestly, like the one in my answer), because the gauge-invariance of the Schrodinger equation guarantees that the energy is the same in any gauge. Thus it is fine to calculate the energy in a gauge where the proton's potential is $\phi=e/r$.
A dissenting commenter seems not to understand that the Schrodinger equation is gauge-invariant. He has argued (in the "Energy/mass of Quantum Vacuum" thread) that rephasing the wavefunction as $\psi \rightarrow e^{i \alpha t} \psi$, where $\alpha$ is constant, shifts the spectrum by $\alpha$. He seems to have forgotten that rephasing the wavefunction is only one part of a gauge transformation. The other part is changing the electromagnetic potentials. When one does both, the energy does not change. This is the entire point of having gauge fields... they are there to make the equations gauge-invariant. His contention in the other thread that "The energy of the ground state is not an observable quantity", apparently because he thinks that it is gauge-dependent, is false.
One potential confusion is that, although the Schrodinger equation for a charge in an electromagnetic field is gauge-invariant, the Hamiltonian is not, in general. But this is not inconsistent with the energy being gauge-invariant, because the Hamiltonian is not always the energy operator. The problem is that a gauge transformation can turn a non-time-dependent Hamiltonian into a time-dependent one, in which case it is no longer true that $\hat{H}\psi=E\psi$.
