Are we overcounting the interaction energy in the classical EM field Lagrangian? The classical electromagnetic field Lagrangian density can be written as
$$-\frac{1}{4} F^{\mu \nu} F_{\mu \nu} - \frac{1}{c} j^\mu A_{\mu} = \frac{1}{2}\vec{E}(\vec{x})^{2} -\frac{1}{2}\vec{B}(\vec{x})^{2} - \rho (\vec{x}) \phi (\vec{x}) + \vec{\jmath}(\vec{x})\cdot\vec{A}(\vec{x}).$$
However, the electric field is subject to the constraint $\vec{\nabla} \cdot \vec{E} = \rho$, and we can write $\vec{E}= - \partial_t \vec{A} - \vec{\nabla}\phi$, essentially splitting it into its transverse (divergence-free) and longitudinal (rotation-free) components. The integral over the field energy of the electric field then yields (because the fields are orthogonal in each Fourier mode):
\begin{align}
\int d^3x \, \left( \vec{E}_{t}\!{}^{2} + \vec{E}_{l}{}^{2} \right) 
&= \int d^3x \, \left[ \vec{E}_{t}\!{}^{2} + (\vec{\nabla} \phi )^2 \right] \\
&= \int d^3 x \, \left(\vec{E}_{t}\!{}^{2}- \phi \vec{\nabla}{}^{2} \phi \right) = \int d^3 x \, \left(\vec{E}_{t}\!{}^{2} - \phi \rho \right).
\end{align}
I have used one of Green’s identities here, and afterwards applied Gauss' law (in the Coulomb gauge).
So the field energy contains (after one applies the constraint (Gauss' law, which is one of the Euler-Lagrange equations)) already the interaction with a charge distribution. What's going on here? Does that mean the other term is counted too much, and is not necessary? Or should I simply not worry, because the Euler-Lagrange equations do follow from this Lagrangian, and can't be applied beforehand? If so, are there any implications connected with this?
If one does the Legendre transform, and writes down the Hamiltonian, one still has 2 terms with the interaction energy, after applying the equations of motion.
 A: The short answer is no, there is no over counting in the standard expressions (of course).
In electrostatics, we can set $\vec{B}=\vec{A}=\vec{J}=0$. Then the action simplifies to
\begin{equation}
S =  \int {\rm d}^4 x \left[\frac{1}{2} \vec{E}^2 - \rho \phi\right] = \int {\rm d}^4 x \left[\frac{1}{2} \left(\nabla \phi\right)^2 - \rho \phi\right]
\end{equation}
Note that the electric field appears with a factor of $\frac{1}{2}$, which was missing from your expression.
As you said, you can use a combination of integration by parts and Gauss's law (which here becomes the Laplace equation, $\nabla^2 \phi=-\rho$), to rewrite the first term, yielding
\begin{eqnarray}
S &=& \int {\rm d}^4 x \left[-\frac{1}{2}\phi \nabla^2 \phi - \rho \phi\right] \\
&=& \int {\rm d}^4 x \left[\frac{1}{2}\rho \phi - \rho \phi\right]\\
&=& \int {\rm d}^4 x \left[-\frac{1}{2} \rho \phi \right] \\
&=& \int {\rm d}t \int {\rm d}^3 x \left[-\frac{1}{2} \rho \phi \right]\\
&=& - E T
\end{eqnarray}
where $T \equiv \int {\rm d} t$ and
\begin{equation}
E = \frac{1}{2}\int {\rm d}^3 x \rho \phi
\end{equation}
is a standard representation of the energy in electrostatics.
So, there is a nice relationship between the action, and the electrostatic energy, in the static case.
Of course, we wouldn't actually expect the action to equal the energy. For example, the action for a point particle in a potential is $S=\int {\rm dt} \left(\frac{1}{2} m \dot{x}^2 - V(x)\right)$, while the energy for a point particle in a potential is $E=\frac{1}{2} m \dot{x}^2 + V(x)$. Without a potential (ie, if $V=0$), then for a point particle, we would find $S=ET$, which is the same as what we found in electrostatics, up to a sign. To properly compute the energy, you should perform a Legendre transform on the electromagnetic Lagrangian to derive the Hamiltonian. This is somewhat subtle due to the existence of a first class constraint associated with electromagnetic gauge invariance. The upshot is that the minus sign relating the action and the energy we found here, is due to the fact that the electric potential $\phi$ is associated with the $0$ component of $A_\mu$, which has the "wrong sign" kinetic energy.
