I claim that $M(n)=n^{\Theta(n)}$.
The upper bound is easy: we can assume without loss of generality that an $n$-line program only uses variables $a_0,\dots,a_{n-1}$, hence there are only $n^2+2n$ possible instructions, and $(n^2+2n)^n$ programs. Thus, one of the numbers $0,\dots,(n^2+2n)^n$ cannot occur as the exit value of such a program, i.e.,
$$M(n)<(n^2+2n)^n=O(n^{2n}).$$
Now, for the lower bound. Note that $M(n)=n^{\Omega(n)}=2^{\Omega(n\log n)}$ can be equivalently restated as saying that every $m$-bit integer can be computed as the exit value of a program with $O(m/\log m)$ lines, so this is what we need to show.
First, it is easy to see that every such number can be computed by a program with $5m$ or so lines (leading to the bound $M(n)\ge2^{n/5}$): the 4-line program snippet
\begin{align*}
\text{loop:}&a_0{-}-\\&a_1{+}+\\&a_1{+}+\\&a_0?\text{ loop}
\end{align*}
computes $a_1:=2a_0$ and clears $a_0$ (provided initially $a_0>0=a_1$), and with one more increment, we can also do $a_1:=2a_0+1$ in 5 lines. By chaining $m$ such snippets (alternating between $a_0$ and $a_1$), we can produce any $m$-bit number.
For $O(m/\log m)$ lines, we have to be more sophisticated. Assume w.l.o.g. that $m=2^k$ is a power of $2$. Then an $m$-bit number can be identified with the truth table of a Boolean function $f\colon\{0,1\}^k\to\{0,1\}$.
Such Boolean functions can be computed by Boolean circuits. It will be convenient to represent here circuits as straight-line programs: a circuit $C$ of size $s>k$ computes $s$ Boolean values $a_1,\dots,a_s$ where $a_1,\dots,a_k$ are initialized to the $k$ input bits of $f$, and for each $k<i\le s$, we have an instruction of one of the forms
\begin{align*}
a_i&:=a_j\land a_k,\\
a_i&:=a_j\lor a_k,\\
a_i&:=\neg a_j,
\end{align*}
where $j,k<i$. The final value of $a_s$ is the output of $f$.
We can compute the number whose binary representation is the truth table of $f$ by an SL-program with the following structure:
\begin{align*}
\text{loop:}&\textit{/* assume $a_1,\dots,a_k$ hold an intended input of $f$: */}\\
&\text{simulate the computation of $C$}\\
&\text{double $a_0$}\\
&\text{if $a_s>0$, $a_0{+}+$}\\
&\textit{/* increment $a_1,\dots,a_k$ as a $k$-bit binary integer: */}\\
&\text{if $a_1=0$: $a_1:=1$, go to loop}\\
&a_1:=0\\
&\text{if $a_2=0$: $a_2:=1$, go to loop}\\
&a_2:=0\\
&\dots\\
&\text{if $a_k=0$: $a_k:=1$, go to loop}\\
&\textit{/* halt */}
\end{align*}
It is easy to see that each instruction of $C$, as well as each of the remaining lines of the pseudocode above, can be implemented with a constant number of SL instructions, hence the total length of the program is $O(s)$.
Now, the crucial point is that by a nontrivial result in circuit complexity going back to Shannon, every Boolean function in $k$ variables can be computed by a circuit of size $s=(1+o(1))2^k/k$ (and this bound is tight for the vast majority of Boolean functions). Thus, every number with $m=2^k$ bits can be output by an SL-program of length $O(2^k/k)=O(m/\log m)$, as claimed.
EDIT: One way of proving the weaker (but sufficient above) bound that every $f\colon\{0,1\}^k\to\{0,1\}$ is computable by a circuit of size $O(2^k/k)$ is as follows. First, by recursively expanding
$$\tag{$*$}f(x_1,\dots,x_k)=(x_k\land f(x_1,\dots,x_{k-1},1))\lor(\neg x_k\land f(x_1,\dots,x_{k-1},0)),$$
we see that $f$ has a circuit of size $2k+3(2^k-1)$ or so. We can shorten it by observing that this circuit has many redundancies: there are nodes computing the function $g(x_1,\dots,x_d)=f(x_1,\dots,x_d,a_{d+1},\dots,a_k)$ for each $d<k$ and $a_{d+1},\dots,a_k\in\{0,1\}$, and many of these functions actually coincide. We can exploit this by precomputing the values of all $2^{2^d}$ Boolean functions in variables $x_1,\dots,x_d$. This can be done by a circuit of size $2^{2^d}$: take the concatenation of arbitrary circuits computing all these functions, and remove redundant nodes computing a function that is also computed by another node earlier in the circuit. After this reduction, no two nodes in the circuit compute the same function, hence there are only $2^{2^d}$ nodes.
Now, if we apply the expansion $(*)$ only until we reach functions in $d$ variables, and use the precomputed values for these, we obtain a circuit for $f$ of size
$$s=2(k-d)+3(2^{k-d}-1)+2^{2^d}.$$
Taking $d=\lfloor\log_2(k-2\log_2k)\rfloor$, we have
\begin{align*}
2^{k-d}&\le\frac{2^{k+1}}{k-2\log_2k}=(2+o(1))2^k/k,\\
2^{2^d}&\le2^k/k^2=o(2^k/k),
\end{align*}
hence $s\le(6+o(1))2^k/k$.
