This is frustrating because there is a lot of study of this sort of question in the combinatorial commutative algebra literature, and the exact example you discuss is a favorite example in this field, but I can't find an answer to your question. So, here are some pointers to the relevant background.

Fix $n$. Let $R$ be the semigroup ring corresponding to the semigroup of $n \times n$ nonnegative integer matrices all of whose row and column sums are equal. So the Hilbert series of $R$ is
$$h(x) := \sum_k |M(n,k)| x^k$$
By the Erhart theorem, $|M(n,k)|$ is a polynomial in $k$, so we can write
$$h(x) = \frac{\delta(x)}{(1-x)^{(n-1)^2+1}}$$
for some polynomial $\delta(x)$ of degree $\leq (n-1)^2$. Richard Stanley pioneered the study of the relation between commutative algebra properties of $R$ and combinatorial properties of $\delta$. In particular, the ring $R$ is Gorenstein (the canonical module is generated by the all ones matrix) and this implies that $\delta$ is palindromic with positive coefficients.
The standard reference for this material is Stanley's book "Combinatorics and Commutative Algebra"; it does not appear to be legally available online.

It might be worth pausing for an example: According to this webpage,
$$|M(3,k)| = 1 + \frac{9 k}{4} + \frac{15 k^2}{8} + \frac{3 k^3}{4} + \frac{k^4}{8}$$
and we can compute
$$h(x) = \frac{1+x+x^2}{(1-x)^5}.$$
The polynomial $\delta$ is $1+x+x^2$.

Now, we have the following implications:

(1) $\delta(x)$ has all real roots $\implies$

(2) Writing $\delta(x) = \sum \delta_k x^k$, we have $\delta_k^2 \geq \delta_{k-1} \delta_{k+1}$ $\implies$

(3) We have $M(n,k)^2 \geq M(n,k-1) M(n,k+1)$ $\implies$

(4) $M(n,k) M(n,k+n-1) \geq M(n,k-1) M(n,k+n)$, which is the relation you want.

The implication $(3) \implies (4)$ is elementary; the others are discussed in Stanley's superb survey "Log-concave and Unimodal sequences in Algebra, Combinatorics, and Geometry".

In that Survey, Stanley made Conjecture 4, that the Hilbert series of any Cohen-Macaulay domain should obey (2). $R$ is a Cohen-Macualay domain (any normal semi-group ring is, by a result of Hochster), but Conjecture 4 turned out to be false. According to "Log-concave and Unimodal sequences in Algebra, Combinatorics, and Geometry: an update" by Brenti (scroll down), a counter example can be found in Niesi and Robbiano; I haven't checked this reference. Stanley and Brenti both suggest that the conjecture is more plausible for Gorenstein rings, and a quick skim through the Mathscinet papers which cite them suggest that no Gorenstein counterexample is known.

I went over to Dennis Pixton's webpage which tabulates the Erhart polynomials for this exact problem. The largest value he lists is $n=9$. I computed the corresponding $\delta$ (coefficients available on request) and found that it violated (1) but obeyed (2).

Specifically, for $n=9$, the polynomial $\delta$ has degree $56$. Of the roots, $52$ were real and clearly isolated. There are also four roots at $(-170629.9 \pm 70111.4 i)^{\pm 1}$. (This all assumes you trust *Mathematica*'s numerical algorithms.)

Finally, you should probably know a few keywords: The polytope of $n \times n$ matrices with nonnegative entries and row and column sums equal to $1$ is the "Birkhoff polytope". The more general case where you just fix $(a_1, a_2, \ldots, a_n)$ and $(b_1, b_2, \ldots, b_n)$ with $\sum a_i = \sum b_i$ and ask for row sums $a_i$ and column sums $b_i$ is the "transportation polytope".