The missing $1/c$ in your first expression is simply a consequence of the units used. The second expression is in Gaussian units while the first one is in either SI units or in natural units. In the latter system of units (natural units) certain constants like $\hbar$ and $c$ have a numerical value of 1, so they can be left out of the equations.$^1$ This is common practice in physics and it doesn't change anything about the dimensional analysis of the problem, as long as you keep in mind that you're working with those natural units.
The same goes for any other system as well. Every system of units $A$ is consistent with any other system of units $B$ as long as you yourself are consistent in their usage and correctly transform everything between $A$ and $B$ when desired.
So there is no fundamental difference between a dimensional analysis in SI, Gaussian or natural units, as long as you keep in mind what units you're working with. The units themselves will (obviously) vary between systems, but dimensional analysis in one system will be entirely consistent with dimensional analysis in another.$^2$
$^1$ Note that this is not the case for SI units. As is rather well-known, the numerical value of $c$ in SI units is about $3\times10^8$ $(\mathrm{m/s})$. The reason for the absence of $1/c$ in SI units is a conventional difference. Wikipedia has a comparison between Gaussian and SI units explaining the major differences here.
$^2$ Perhaps one important note concerning Gaussian and SI units here is that due to the different conventions, it can be more difficult to transform between them. E.g. making an equation dimensionless in SI units, might yield a non-dimensionless equation when transformed into Gaussian units.
One example is when we consider Gauss's law in Gaussian units divided by the free charge density: $(1/\rho)\vec{\nabla}\cdot\vec{E} = 4\pi$. The quantity on the left-hand side is dimensionless in Gaussian units, but not in SI units, where it is $(1/\rho)\vec{\nabla}\cdot\vec{E} = 1/\epsilon_0$. So you have to watch out for that when transforming your equations. Dimensional analysis may therefore also yield seemingly different results in SI or Gaussian units, but there is no problem if you remember the conventional differences and, again, stay consistent.
I will try to slightly elaborate on @VladimirKalitvianski answer.
From Maxwell's equations, we can derive that the following combination of gauge transformations on $\mathbf{A}$ and $\Phi$ leave both $\mathbf{B}$ and $\mathbf{E}$ invariant:
\begin{align}
& \mathbf{A}'=\mathbf{A}-\mathbf{\nabla} \alpha \\& \Phi'=\Phi+\frac{\partial \alpha}{\partial t}
\end{align}
where $\alpha=\alpha(\mathbf{x},t)$. This means that all the field configurations of $\mathbf{B}$ and $\mathbf{E}$ related by a gauge transformation are physically equivalent. Note that this has nothing to do with the Hamiltonian operator in QM.
Now in QM, we know that a wave function can always be multiplied by a phase factor:
$$ \psi'=e^{-iq\alpha}\psi, $$
where $\alpha \neq \alpha(\mathbf{x},t)$, because the probability of finding the particle at a particular position is unaffected by the above transformation, and also the Schrodinger equation and the probability current are unaffected by the above transformation. If we now demand that the above also holds for when $\alpha = \alpha(\mathbf{x},t)$ (i.e. a gauge transformation), then the Schrodinger equation must be made gauge invariant:
\begin{equation}
i\frac{\partial\psi}{\partial t}=-\frac{1}{2m}(\mathbf{\nabla}-iq\mathbf{A})^2\psi+(V+q\Phi)\psi
\end{equation}
such that the Schrodinger equation is invariant under the simultaneous gauge transformations:
\begin{align}
&\mathbf{A}'=\mathbf{A}-\mathbf{\nabla} \alpha \\&
\Phi'=\Phi+\frac{\partial \alpha}{\partial t} \\&
\psi'=e^{-iq\alpha}\psi \tag{1}
\end{align}
Note that we can say that we have adjusted the "normal" Hamiltonian by replacing the ordinary (partial) derivatives by:
\begin{equation}
\begin{array}{cc}
\displaystyle \mathbf{\nabla} \rightarrow \mathbf{D}\equiv \mathbf{\nabla} - i q \mathbf{A} \; ,& \displaystyle \frac{\partial}{\partial t} \rightarrow D^0 \equiv\frac{\partial}{\partial t} + iq
\end{array}
\end{equation}
To sum up, by demanding that our theory is invariant under the gauge transformation expressed by equation $(1)$, we are forced to change the Hamiltonian operator as we have done above. However, by doing this, the new Hamiltonian describes a particle interacting with the potentials $\mathbf{A}$ and $\Phi$. If you're not convinced by this argument, I strongly recommend you to read up on the Aharonov–Bohm effect (http://en.wikipedia.org/wiki/Aharonov%E2%80%93Bohm_effect).
Furthermore, note that we require that a gauge transformation does not affect any observables. This means that we must demand that the probability current is also unaffected. You can show (although it is quite tedious) that the current is made gauge invariant by making the replacement: $\displaystyle \mathbf{\nabla} \rightarrow \mathbf{D}$.
Best Answer
For a particle-field system the only way to define a gauge invariant energy is to consider the energy carried by the field as well, in the form of the energy momentum tensor $T^{\mu\nu}$ in the presence of charges. $T^{\mu\nu}$ is a manifestly gauge invariant quantity. To derive this use Noether's theorem and the maxwell equations in the presence of charges. This gives a conserved energy for the system
So the answer to this would be no, it does not correspond to the total energy as you must also consider the field energy to attain a conserved quantity. Intuitively, if we choose some gauge where $\phi = 0$, then the change in the energy of particle is compensated for by a change in the energy carried by the fields. Thus resulting in the same total energy for the system.