This is a partial answer for the special case when all $a_i$ are distinct.
I will work on complex(or algebraic closure of $\mathbb{Q}$) for convenience.
Writing $M-\lambda I = D_{\lambda}+E$ where $D_{\lambda}$ is the diagonal matrix with diagonal entries $a_0-\lambda+1, \cdots , a_n-\lambda+1$.
With this expression, it is easy to calculate the determinant, which will give the characteristic polynomial of $M$.

If $f(\lambda)=(a_0-\lambda+1)\cdots (a_n-\lambda+1)$, then we have
$$
\textrm{det}(M-\lambda I) = f(\lambda)+f'(\lambda).$$
Considering the identity
$(f(t)e^t)'=(f(t)+f'(t))e^t$,
we can find the roots of characteristic polynomial by looking at the critical points of $f(t)e^t$.
Also, by Mean Value theorem, we know that the critical points of $f(t)e^t$ are all distinct.

Therefore, if we let $\lambda_0, \cdots, \lambda_n$ be the critical points of $f(t)e^t$, then we have the following Smith Normal form of the matrix $M-\lambda I$ over the polynomial ring $\mathbb{C}[\lambda]$, (or $\overline{\mathbb{Q}}[\lambda]$):
$$
\textrm{Diag}(1,\cdots, 1, (\lambda-\lambda_0) \cdots (\lambda-\lambda_n)).$$
Hence, we obtain the Smith Normal form of $M$ in this case:
$$\textrm{Diag}(1,\cdots, 1,\lambda_0 \cdots \lambda_n).$$
There is a natural way of bringing this down to $\mathbb{Z}$, then we have to deal with the irreducible factors of $f(t)+f'(t)$.
This is indeed
$$\textrm{Diag}(1,\cdots, 1, f(0)+f'(0)).$$
Added) This method also works for the case below:

The cardinality of $\{ j: a_i = a_j \}$ is at most $2$, for any $i=0,\cdots, n$.

Added2) This gives the SNF of $M$ over $\mathbb{Q}$. Over $\mathbb{Z}$ will be certainly more difficult.