On determinants formed by binomial coefficients Let $q$ be a number. Let us consider the $q^2-1$-th line of the Pascal triangle (i.e. numbers ${{q^2-1} \choose i}$, $i=0,1,...q^2-1$). We have $q^2$ numbers. 
Let us form naively a $q \times q$ matrix from them. Example for $q=2$:
\[\begin{pmatrix} 1 & 3 \\
3 & 1 \end{pmatrix}\]
And let us take its determinant: $-8$ if $q=2$. 
How to prove that if $q$ is prime then $q$ enters in this determinant with 
$(q+4)(q-1)/2 $-th exponent?
Example: $q=2$, $(q+4)(q-1)/2 =3$. Really, $8=2^3$. 
Case $q=3$: the binomial coefficients are $1,8,28,56,70,56,28,8,1$.
The matrix: 
\[\begin{pmatrix} 1 & 8 & 28 \\
56 & 70 & 56 \\
28 & 8 & 1 \end{pmatrix}\]
The determinant is $2 \cdot 3^7 \cdot 7$. Really, $(q+4)(q-1)/2 =7$. 
 A: Looking at Krattenthaler's famous "determinant calculus" survey  (also, suggested by Per Alexandersson) we find Theorem 26 in there, which answers a more general question, and should (most likely) yield the claim in the OP.

Theorem Let $q$ be a nonnegative integer, and let $L_1,\ldots,L_q$ and $A, B$ be indeterminates. Then,
  \begin{equation*}
 \det\left( \binom{BL_i+A}{L_i+j}\right) = \frac{\prod_{i < j}(L_i-L_j)}{\prod_{i=1}^q(L_i+q)!}\prod_{i=1}^q\frac{(BL_i+A)!}{((B-1)L_i+A-1)!}\prod_{i=1}^q(A-Bi+1)_{i-1},
\end{equation*}
  where $(a)_k$ is the usual Pochhammer symbol.

The matrix in the OP is given by
\begin{equation*}
 M_{ij} = \binom{q^2-1}{q(i-1)+j-1},
\end{equation*}
thus we can select $A=q^2-1$, $B=0$, $L_i = q(i-1)-1$, and apply the above result, and hopefully things fall into place.
A: For the specific determinant in the question there is a slightly more direct proof than Lemma 3 in Krattenthaler. We can evaluate it using the Jacobi-Trudi identity. Let's denote  $\lambda _i=(q-i+1)(q-1)$. The determinant in question, $\det\left(\binom{q^2-1}{qi-q+j}\right)$, is equal (up to sign) to the determinant $\det\left(e_{\lambda_i-i+j}\right)$ when you put all variables equal to $1$. But this is equal to the Schur polynomial of the conjugate partition $\lambda^{\star}=(q^{q-1}(q-1)^{q-1}\cdots 1^{q-1}0^{q-1})$ evaluated at $1,1,\dots,1$. (I added a bunch of zeros to the partition so that the number of parts matches the number of variables.)
There is a simple product formula for Schur polynomials evaluated at $1,1,\dots,1$, and in this case it gives
$$s_{\lambda^{\star}}(1,1,\dots,1)=\prod_{i<j}\frac{\lambda^{\star}_i-\lambda^{\star}_j+j-i}{j-i}.$$
The denominator has the form $1!2!\cdots (q^2-2)!$ so the highest power of $q$ dividing it is $q^{\frac{(q-1)(q^2-2)}{2}}$. For the numerator, notice that $\lambda^{\star}_i-\lambda^{\star}_j+j-i$ is divisible by $q$ iff $j-i$ is divisible by $q-1$. This means that exactly $(q-1)\binom{q+1}{2}$ factors in the numerator are divisible by $q$. Out of these the pairs $(i,j)=(r,q^2-q+r)$ for $r\in\{1,2,\dots,q-1\}$ are divisible by $q^2$. So the power of $q$ dividing your determinant is exactly
$$(q-1)\binom{q+1}{2}+(q-1)-\frac{(q-1)(q^2-2)}{2}=\frac{(q-1)(q+4)}{2}$$
as desired.
