If the vector $\mathbf X$ has a multivariate Gaussian distribution $\mathcal N(\mathbf \mu, \mathbf \Sigma)$ and the random variable $\mathbf Y$ is defined as $\mathbf Y = \mathbf M \mathbf X$, where $\mathbf M$ is a constant matrix, then $\mathbf Y$ has a multivariate Gaussian distribution $\mathcal N(\mathbf M \mathbf \mu, \mathbf M \mathbf \Sigma \mathbf M^T)$. (See this Wikipedia page.)
In your situation, $\mathbf Y = (X_1 - X_2, \dots, X_1 - X_M)$ is related to $\mathbf X = (X_1, X_2, \dots, X_M)$ by the equation $\mathbf Y = \mathbf M \mathbf X$, where
$$ \mathbf M = \begin{bmatrix} 1 & - 1 & 0 & \dots & 0 \\ 1& 0 & - 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots& \vdots\\1 & 0 & 0 & \dots & -1\end{bmatrix}$$
Judging from the answer you linked, I suspect you want
$$ \mathbf \mu = \begin{bmatrix} \mu \\ \mu \\ \vdots \\ \mu \end{bmatrix}, \ \ \ \ \ \mathbf \Sigma = \begin{bmatrix} 1 & \rho & \dots & \rho \\ \rho & 1 & \dots & \rho\\ \vdots & \vdots & \ddots & \vdots \\ \rho & \rho & \dots & 1\end{bmatrix}$$
So by the result above, $\mathbf Y$ has a multivariate normal distribution, with mean
$$ \mathbf M \mathbf \mu = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0\end{bmatrix}$$
and covariance
$$ \mathbf M \mathbf \Sigma \mathbf M^T = \begin{bmatrix} 2(1- \rho) & 1 - \rho & \dots & 1 - \rho \\ 1 - \rho & 2(1 - \rho) & \dots & 1 - \rho \\ \vdots & \vdots & \ddots & \vdots \\ 1 - \rho & 1 - \rho & \dots & 2(1 - \rho)\end{bmatrix} $$
[And by the way, in the answer you linked for the $M = 2$ case, $\Sigma$ is the matrix $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$; it is not $I \times \rho$. The constant $\rho$ can be interpreted as the correlation between $X_1$ and $X_2$.]