How can the first Maxwell equation be valid in non-static cases? I am thinking in the framework of Classical Gravity, where the speed of the interaction is infinite. Now it is also known that there is a correspondence between Classical gravity and electrostatics, in fact the analogous of the first Maxwell equation can be derived assuming Newton law for a pointlike mass and the superposition principle, giving us:
$$
\nabla \cdot {\bf{g}} \, = \, - 4\pi G \rho_M({\bf{x}}) \tag{1}
$$
Where ${\bf{g}}$ is the gravitational field and $\rho_M$ is the matter density.
This is analogous to the electrostatic (but valid in General) Maxwell equation, which is, in gaussian units:
$$
\nabla \cdot {\bf{E}} \, = \, - 4\pi \rho_E({\bf{x}}) \tag{2}
$$
I really don't understand how can the second equation holds when charges are in motion, let me try to explain.
The fields ${\bf{g}}$ and ${\bf{E}}$ satisfy the same equation, and at the same time can be made arbitrarily (?) different because in Newtonian mechanics the speed of the interaction is infinite and therefore the ${\bf{g}}$ of a pointlike mass is always radial, while the Electric field of a moving charge can be as messy as you want it, depending how you move it but somehow the ${\bf{E}}$ it still satisfies the very same relation the ${\bf{g}}$ field satisfies.
I mean how can The Electric field rearrange itself instant by instant such that it satisfies the same relation of a completely different(?) gravitational field (static electric field) ?
 A: When you're looking at the differential form of Maxwell's equations, there's no need to worry about propagation speed. The quantity $\rho(\mathbf{x})$ on the right-hand side only depends on what's at the point $\mathbf{x}$, and so does the left-hand side. The equation is manifestly local, so information doesn't have to propagate anywhere; it's all at $\mathbf{x}$.
So there's no problem with the differential forms, but you might worry about the integral forms, which look "nonlocal". Really, there's no problem with them either, because they're completely mathematically equivalent to the differential forms, but it's worth thinking a bit about why they work. For example, Gauss's law in integral form is 
$$\oint \mathbf{E} \cdot d \mathbf{S} = \frac{Q_{\text{in}}}{\epsilon_0}.$$
You might think that if the surface used on the left-hand side is big, then it can't react instantly to changes in $Q_{\text{in}}$. But there's no problem because charge is locally conserved, which you can prove using the rest of Maxwell's equations. For $Q_{\text{in}}$ to change, charge can't just appear out of nowhere; it has to cross the surface, and that's how the flux through it "knows" to change. 
A slightly more puzzling question is how Faraday's law can work,
$$\oint \mathbf{E} \cdot d \mathbf{s} = - \frac{d \Phi_B}{dt}$$
when there is no corresponding "conservation law" for magnetic flux. I address that here. Gauss's law for magnetism and Ampere's law can be addressed similarly. 
A: I think you might have made a really good observation. The finite speed of light is "hidden" in the 4 Maxwell Equations. That means you need both the divergence and curl. 
Jefimenko's Equations explain how the fields might satisfy certain relations instantaneously. The fields don't actually induce time variation in each other, they are generated by the same time varying sources and those dictate the fields' time variation which happens to satisfy patterns with respect to each other. I don't think there's an equivalent for gravity. 
One problem with the gravitational form of Gauss' Law is that the gravitational field/space time curvature is not completely determined by mass density. It also relies on energy and stress/pressure. Thisr introduces a time varying component. So to some degree, ignoring time variation is built in. I'm not sure how to recover that without a full on GR treatment. 
I stand corrected there was some effort to Maxwellize Newton before Special Relativity: Gravitomagnetism
A: Part of the question is : how can it be that $\vec \nabla \cdot \vec E= -4\pi \rho$ is obeyed when charges are moving. Indeed the equation appears to indicate infinite speed of propagation. However when charges are moving E has a rotation. If this is taken into account it turns out that E propagates at speed $c$. In the Lorentz gauge $(A^0,\vec A)$, $\vec E$ and $\vec B$, obey inhomogeneous wave equations with phase speed $c$, as is conveniently shown when the covariant notation is used. I leave this as an exercise to the reader.
