Distribution of determinants of $n\times n$ matrices with entries in $\{0,1,\ldots,q-1\}$

Consider the set $M(n,q)$ of $n\times n$ matrices with entries in $\{0,1,\ldots,q-1\}$, where $q$ is a prime power. What can be said about the distribution of the determinant of matrices in $M(n,q)$? (A 'heuristic' statement of the problem: taking $\{0,1,\ldots,q-1\}$ as a basis for $F=\mathbb{Z}_q$, what do the determinants of matrices over $F$ look like if you don't mod out $q$?)

Obviously $|M(n,q)| = q^{n^2}$. Since $|GL_n(\mathbb{F}_q)| = \prod_{k=0}^{n-1} q^n-q^k$, in $\mathbb{F}_q$ we get a clean answer for how many are divisible by $q$: the values are equally distributed (modulo $q$, there are $\frac{1}{q-1}\prod_{k=0}^{n-1} q^n-q^k$ matrices with determinant $j$, $1\le j\le q-1$). But if we do not look mod $q$, as it were, the question becomes substantially more difficult; to be frank, I'm not sure where to start or if there are any clear patterns. Information about the limiting behavior or any upper bounds on the magnitude of the determinant would be welcome as well.

I computed the distributions for several values of $n=2,3$ and $2\le q\le 5$; the plot labels are of the form $\{n,q\}$.

As expected, determinant zero is the most common option and a determinant of $a$ is just as likely as a determinant of $-a$. Past that, I admit I'm a little out of my league, but it seems like an interesting problem.

In the continuous limit of $n$ fixed, $q \to \infty$, $\frac{1}{q^n}\log |\det(M)|$ is asymptotically normal as $n \to \infty$. See Terry Tao's comments on this MO thread. The linked paper of Nguyen--Vu has a nicely readable intro, see particularly around equations (1.6)-(1.7). The intuition is roughly that the determinant is going to be the (signed) hypervolume, which can be computed as an iterated multidimensional "base times height". Taking the logarithm and fuzzing your eyes, this looks like a sum of i.i.d. random variables. The details are involved, of course, and I have not attempted to digest them. Someone with expertise in this approach may be able to "discretize" it quickly. It's likely more appropriate as an MO question than an MSE question.

Edit: Now that I look at it, Nguyen--Vu's actual main Theorem 1.1 covers this discrete case too (even though the linked MO thread was just after the continuous case), and more generally any distribution with exponentially decaying tails. So, we get that for fixed $q$, $\log |\det(M)|$ is asymptotically normal as $n \to \infty$, with explicit convergence rates.