Why doesn't Gauss's law for gravity apply for an unbounded, continuous, and homogeneous mass? Consider Gauss's law for gravity, in its differential form:
$$\vec{\nabla}\cdot \vec{g}=-4\pi G\rho,$$ or in its integral form:
$$\iint\vec{g}\cdot d\vec{A}=-4\pi G M.$$
This law intuitively makes sense. If no matter is present no net flux is going through a surface enclosing a volume if a mass is present in this volume.
But what if we consider an unbounded, continuous, and homogeneous distribution of mass? The value of $\vec{g}$ is zero everywhere in the distribution. So the surface integral of $\vec{g}$ over every closed surface is zero. This clearly contradicts Gauss's law, which says that this integral has the value $-4\pi GM$, $M$ being the total mass ($\rho dV$) inside the volume.
why is that so? Is there a mathematical reason (though I can't see why this should be the case) or is the existence of the proposed mass distribution just not realistic? Or what?
I make an edit because I think I made a mistake. Concerning the interpretation of Gauss's law. In evaluating the surface integral (so I think), one has to take the value of the $\vec{g}$ field at the surface. One must only consider the $\vec{g}$ field produced by the mass inside (and on) the surface. Which is perfectly well defined! When one has a distribution of electric charges one can encapsulate some of them and calculate the Gass integral. The integrals of the other charges cancel (they are outside the shell). As long as the distribution is finite. If the distribution is unbounded, you can't say anymore that the fields produced cancel. Simply because the field they produce on the shell enclosing a bound part of the distribution is unbounded (infinite), and in this case, the integral is not well defined. Unbounded values (infinities) just don't cancel when summed in every direction.
Somehow I think this question is connected with this question that was asked two days ago. The question asks if it's allowed to consider continuous charge distributions in Maxwell's equations. Can we consider continuous mass distributions in classical Newtonian gravity? If we consider discrete mass distributions we can imagine a surface around every mass, or collection of masses. The contribution of the outside masses to the total gravitational flux through the surface is zero. Only the masses inside the surface contribute to the total flux. The integral form of Gauss's law refers to flux (surface integral), while the differential form refers to local (point) values. On a planet in a universe with an unbounded amount of discrete masses, I will still be able to stand, i.e., there is a force of gravity. The flux through the surface surrounding the planet, due to all other masses, will be zero though. Or not?
But what if we make the distribution continuous? I won't be able to stand on the planet anymore (apart from the fact that I can't walk through a continuous mass distribution). There is zero gravity. The contribution to $\vec{g}$, on a point of the enclosing surface (of the planet, which has now become a sphere of continuous mass surrounded by an unbounded amount of continuous mass), will be the sum of two unbounded contributions: that from the unbounded mass residing in the space on one side of the tangent plane to the point, and one on the other side. They will cancel $\vec{g}$, for each point on the surface) caused by the mass inside the surface. Although the two contributions are infinite, the difference will be finite. I vaguely see a connection with an affine space (and with renormalization in quantum field theory: here infinite masses are rendered finite, but one can also keep the masses infinite and just look at mass difference; ohooooh, what I've written?).
So Gauss's law doesn't apply for continuous, unbound masses (ohoooh, what I've written? If we make the continuous distribution discrete (an unbounded number of separate continuous mass distributions), then we can walk on every mass in this distribution (if they are big enough). And we can place an enclosing surface around each mass and see that the force of gravity is non-zero at the surface of each mass.
Ohooooh! I'm sure I've overlooked something (mathematical rigor?).
The integral form of Gauss's law is (for me) much easier to digest than the differential form (which is the one that allegedly caused a problem in the linked question). For example in empty space, the integral form is obvious, but what about the differential form? Indeed, $\rho$ is zero out there, but $\vec{g}$ doesn't have to be.
 A: The problem is that you are proposing infinite forces without realizing it. Consider any point in the infinite homogeneous mass and draw a plane through that point. The total gravitational force on a small lump of matter at the chosen point, from all the material on one side of the plane (or if you prefer, from all the material within a region subtending some given finite solid angle at the point), is infinite. The total gravitational force from the material on the other side of the plane is also infinite. The total force is the sum of these infinite opposing forces. That sum is not zero, but undefined.
From a physical standpoint one may say that the problem has been posed in an unphysical way: one has proposed an impossible physical system, and Newtonian gravity has no prediction to make other than to warn us that the situation is indeed unphysical. From a mathematical standpoint, the conclusion is that arguing from symmetry has to be approached carefully when infinity is in play. An equality that holds for finite quantities may not be mathematically well-defined for infinite quantities.
As soon as you make the distribution finite, the forces also become finite, and now $\bf g$ is well-defined and everything works out correctly.
Added remark about homogeneity
I would like to add that the assumption of perfect homogeneity is also questionable here. In science generally, it is often valid to assume homogeneity as a first approximation. In order to do this one need not think the system is really perfectly homogeneous, but a homogeneous model is useful because it approximates closely enough to what the system is really like that it can be used to gain good insight. However, in an infinite system with long-range effects (e.g. gravitational or electromagnetic) a tiny departure from homogeneity, over a large enough region, could cause large or even infinite forces. And since in fact a system will have departures from homogeneity at some level, in this case one may not assume that a homogeneous model yields a good first-order approximation. It may, or it may not: further investigation is needed.
A: The value of ${\bf g}$ will not be constant in this setting.  Nor will ${\bf g}$ be well defined with only the infinite charge distribution given. Some boundary conditions are needed. These will depend on how you build up the infinite system as a limit of finite systems.  This very issue crops up in attempts to construct a "Newtonian" version of  FRW cosmology. In this attempt it  is necessary to choose a point to be the center of the universe.
A: While some of the existing answers (of which I find Mike Stone's to be the most enlightening) make valid points, there does not seem  to be enough emphasis on the reason why  the application of Gauss's law fails for OP and some of the answerers. And it's not because Gauss's law “does not work” for infinite mass distributions (Gauss's law is a local equation connecting variation of gravitational field near given point with mass density at that point, there is no reason why it should fail for homogeneous density) but because of the failure to find the most general form of gravitational field compatible with the symmetries of the problem (homogeneity and isotropy). And in order to find such a form we must take into account the inherent ambiguity in defining gravitational field for such infinitary mass system. In turn, the reason for such ambiguity is simple: there are no inertial reference frames in the world under consideration.
Note, that gravity acts on everything, including observers, and if there is constant mass density everywhere, then it is impossible to construct (even in principle) a reference body on which no gravitational force acts, so that it could serve as a foundation for an inertial reference frame. When we consider only finite mass distributions, an “observer at infinity“ (understood as a limit for a sequence of observers more and more distant from all the masses) could serve as such a reference, defining a privileged class of inertial reference frames. But for infinitary system of masses this construction would fail, and all possible reference frames must be considered non-inertial.
While all the reference frames we can construct would be non-inertial, we can at least select non-rotating frames where no Coriolis forces occur. Thus we arrive at a class of preferred, non-rotating reference frames for this problem: their Cartesian coordinate systems $(t,x^a)$ are all related by the transformations:
$$
{x^a}' = R^a{}_bx^b+{d^a}'(t) \tag{t},
$$
where $R^a{}_b$ is a constant rotation matrix, and $d^a(t)$ is a translation vector depending arbitrarily on a time $t$. All physical quantities transform in the obvious ways under such transformation. In particular, even after fixing origin and orientation of a reference frame, gravitational acceleration field at a given moment can only be defined up to an overall additive constant, acceleration of non-inertial reference frame. In other words, there is a “gauge transformation”  for gravitational acceleration:
$$
 \mathbf{g}'(\mathbf{r})= \mathbf{g}(\mathbf{r})+\mathbf{a}. \tag{g}
$$
Gravitational acceleration is thus inherently ambiguous, which is not surprising, since if the whole world is falling freely with an acceleration $a$ everywhere there is no experiment that could detect such an acceleration. Note, that while position, velocity and acceleration by themselves are ambiguous, relative position, relative velocity and relative acceleration between pairs of bodies are not.
Not taking into account the existence of such ambiguities in defining gravitational acceleration explains OP's (and some of the answerer's) failure to find gravitational acceleration field compatible with homogeneity and isotropy of the universe,  since the only vector field that remains unchanged under arbitrary rotations and translations is trivial.

The value of $g⃗$  is zero everywhere in the distribution.

This is the error, because for a gauge field to be compatible with space–time symmetries their action does not have to leave it unchanged, but rather this action must be equivalent to some gauge transformation.  It is easy to see, that for a vector  field:
$$
\mathbf{g}(\mathbf{r})=\alpha \, (\mathbf{r} - \mathbf{b} ),\tag{f}
$$
where $\alpha$ is a scalar parameter, spatial rotations and translations are equivalent to a change in a vector constant $\mathbf{b}$ which is just a gauge transformation $(g)$. It is easy to see that $(f)$ is the most general form of field $\mathbf{g}$ (at a given moment of time) compatible with homogeneity and isotropy.
By inserting the field $\mathbf{g}$ in the form $(f)$ into Gauss's law we find
$$\alpha=-\frac{4\pi}{3}G\rho_0.$$
If we do a gauge fixing by selecting a “center of the Universe” and placing it at the origin, the resulting gravitational field would then be:
$$
\mathbf{g}(\mathbf{r})=-\frac{4\pi}{3}G\rho\, \mathbf{r}.
$$
In order to find the evolution with time of all the quantities we must first introduce initial velocities (also compatible with homogeneity and isotropy and  having a similar ambiguity related to $(t)$), integrate the masses motion under the gravitational field and solve the continuity equation:
$$
\frac{\partial \rho }{\partial t} +\mathbf{\nabla}\cdot (\rho\mathbf{v})=0.
$$
The resulting solution would essentially be equivalent to FRW cosmological solution of general relativity with “dust-like“ matter (up to some parameters redefinition).
A: Since the gravitational field is irrotational, $\vec g$ can be derived from a potential:
$$\vec g = -\vec \nabla \Phi$$
The divergence equation then becomes
$$\Delta \Phi = 4\pi G\rho$$
Laplacian-type PDE's are often solved by the (static) "plane wave" eigen functions of the Laplacian, namely
$$\psi_{\vec k}(\vec r)=e^{i\vec k\cdot \vec r}$$
This amounts to representing $\Phi$ and $\rho$ as Fourier integrals:
$$\Phi(\vec r)=\int \Phi(\vec k) e^{i\vec k\cdot \vec r}d^3 \vec k$$
$$\rho(\vec r)=\int \rho(\vec k) e^{i\vec k\cdot \vec r}d^3 \vec k$$
In the PDE this leads to
$$-\int \vec k^2\Phi(\vec k) e^{i\vec k\cdot \vec r}d^3 \vec k=4\pi G\int \rho(\vec k) e^{i\vec k\cdot \vec r}d^3 \vec k$$
Equating coefficients on the LHS and RHS leads to
$$\Phi(\vec k)=-4\pi G\frac{\rho(\vec k)}{\vec k^2}$$
The coefficients $\rho(\vec k)$ can be obtained by Fourier transform, i.e.
$$\rho(\vec k)=N\int \rho(\vec r) e^{-i\vec k\cdot \vec r}d^3 \vec r$$
with the Fourier normalization coefficient $N$ (I can't remember it everytime I use FT, but it's irrelevant for the conclusion anyway). But you have required $\rho(\vec r)=\rho_0=const.$, and so
$$\rho(\vec k)=\rho_0 N\int e^{-i\vec k\cdot \vec r}d^3 \vec r=\rho_0 \delta(\vec k)$$
This is intuitive because the Fourier component for $\vec k=0$ is nothing else but the constant part of the function to be transformed (in this case the density, which has been set constant by definition). Hence, we have for the static gravitational potential the "solution"
$$\Phi(\vec r)=-4\pi G\int \frac{\rho_0 \delta(\vec k)}{\vec k^2} e^{i\vec k\cdot \vec r}d^3 \vec k$$
This is obviously divergent for $\vec k\to 0$, so no solution exists, at least not in the $L^2$-Hilbert space formalism.
Your assumption that $\vec g=0$ is a solution is wrong, which could also be checked much quicker by substituting it into the Gauss law equation directly:
$$\vec \nabla \cdot \vec g=\vec \nabla \cdot \vec 0=0 \not= -4\pi G\rho_0$$
Even a general constant $\vec g$ would not change this, because its derivatives are all zero as well.
A: To derive Gauss's Law for Gravity, you start with Newton's Law of Gravitation, and integrate over the entire region. Then, you take the divergence of both sides of the equation, and rearrange.
The integral, in this case, might be written like this:
$$
\vec{\bf g}=-G\rho\iiint \frac{{\bf r}-{\bf s}}{||{\bf r}-{\bf s}||^3}\ d^3s
$$
The problem, here, is that the triple integral is an improper integral that does not have a distinct, unique value. Instead, its value changes depending on how you perform the integral - it is conditionally convergent.
To see this, consider the 1D equivalent* of the integral for $r=0$:
$$
\int_{-\infty}^\infty \frac{s}{|s|}\ ds
$$
If you evaluate this around $s=0$, you get
$$
\lim_{a\to\infty} \int_{-a}^a \frac{s}{|s|}\ ds = \lim_{a\to\infty} \int_{-a}^0 -1\ ds + \int_0^a\ ds = \lim_{a\to\infty} -a+a=0
$$
However, if you evaluate it around $s=1$, you get
$$
\lim_{a\to\infty} \int_{1-a}^{a+1} \frac{s}{|s|}\ ds = \lim_{a\to\infty} \int_{1-a}^0 -1\ ds + \int_0^{a+1}\ ds = \lim_{a\to\infty} (1-a)+(a+1)=2
$$
As a result, you can get a different value depending on the choice of reference point. The same applies to our triple integral.
To get $\vec{\bf g}=0$, we need to choose the reference point ${\bf s}={\bf r}$ (which is the physically most reasonable choice). However, you can successfully integrate in spherical coordinates using ${\bf s}=0$ as the reference point, if you perform the integration in a particular order... which produces the solution
$$
\vec{\bf g}=-\frac{2\pi^2 G\rho}{r}\hat{\bf r}
$$
This also isn't a solution to the Gauss's Law equation, as the divergence is $-\frac{2\pi^2 G\rho}{r^2}$. And if you change the specific method of integration, it may change the result further.
The derivation of Gauss's Law assumes that the integral is well-behaved. If it is not, then moving the divergence from outside the integral to inside it cannot be done, and the derivation fails.
When working with a finite total mass, however, the resulting integral is absolutely convergent - the choice of reference point is irrelevant, and you get the same result for all values. In this case, you can move the divergence from outside the integral to inside, and everything works nicely.

* This is the 1D equivalent because $\text{div} \frac{x}{|x|}{\bf i} = 2\delta(x)$ in the same way that $\text{div} \frac{\bf r}{|{\bf r}|^3}=4\pi\delta({\bf r})$.
A: 
Is there a mathematical reason?

Yes.  The equation $\nabla \cdot \mathbf g(\mathbf r) = -4\pi G\rho_0$ on the domain $\mathbb R^3$ does not admit a unique solution in the absence of boundary conditions, even if we also add the equation $\nabla \times \mathbf g = 0$.  After all, $\mathbf g_1(\mathbf r) =( -4\pi G\rho_0 x )\hat x$ is a solution to both equations, as is $\mathbf g_2(\mathbf r)=(-4\pi G \rho_0 y) \hat y$.
In order to distinguish them, we need to prescribe boundary conditions. On an infinite domain, this is done by specifying the limit of the vector field as $\mathbf r$ goes to spatial infinity in every direction$^\dagger$.  Such a prescription would take the form of a vector-valued function on the two-sphere, $\mathbf f: S^2 \rightarrow \mathbb R^3$, which associates a "limit" vector with each direction in space.
All of this is fine, but you would now like to impose the symmetry demands of homogeneity and isotropy on $\mathbf g$, and this is not possible.  Homogeneity implies that $\mathbf g$ is constant, but isotropy implies that $\mathbf g\mapsto -\mathbf g$ after a $180^\circ$ rotation about any axis.  The only possible choice compatible with both symmetries is $\mathbf g=0$ everywhere, but this is not compatible with $\nabla \cdot \mathbf g=-4\pi G\rho_0\neq 0$. Therefore, while there are an infinity of solutions to the equation $\nabla \cdot \mathbf g=-4\pi G \rho_0$ with $\rho_0\neq 0$, $\mathbf g=0$ (the only homogeneous and isotropic vector field) is not one of them.

$^\dagger$This limit need not exist - one could imagine a vector field whose magnitude grows without bound with increasing distance from some reference point.  In such cases, it would suffice to choose an arbitrary surface enclosing your region of interest and prescribing the boundary conditions on that surface.
