For A fixed $M$, such finite index subgroups may be categorized as follows. Let $W\in GL(2,\mathbb{Z})$ and $b$ a positive integer dividing the lower left element of the matrix $W^{-1} M W$. Then the columns of the matrix $W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}$ generate a subgroup $H$ of $\mathbb{Z}^2$ (of finite index $b$) for which $Mv\in H$ for any $v\in H$. Conversely, all such subgroups for this $M$ arise this way, up to scaling the subgroup by an integer $a$ as already noted in the question.

To see this, first note that a finite index subgroup of $\mathbb{Z}^2$ may be generated by the columns of a non-singular integer matrix $A$, known as a basis. Any other basis equals $AU$, where $U$ is unimodular, i.e. an element of $GL(2,\mathbb{Z})$.

For unimodular $M$, you want to know for which such $A$ do we have $MA = AV$, for some integer matrix $V$, as such $V$ just makes integer linear combinations of the columns of the basis $A$, which give elements of the finite index subgroup $H$. But determinants then imply $V$ is also unimodular, as it has the same determinant as $M$. So $M$ takes any basis of $H$ to another basis of $H$.

Integer matrices have a Smith Normal Form: there exist unimodular $W,V$ so that $W\begin{bmatrix} a & 0\\0 & ab\end{bmatrix}V = A$, where $a,b$ are integers (here positive, as the determinant is non-zero). As any choice of basis will do, we choose the basis $W\begin{bmatrix}a & 0\\0 & ab\end{bmatrix}$.

So an $H$ may be specified by a unimodular $W$ and integers $a$ and $b$. As $a$ simply scales everything we henceforth assume WLOG that $a=1$.

We now characterize such pairs $W, b$. We have from before the equation
$$ M (W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}) = (W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}) V$$
so
$$ (W^{-1} M W) \begin{bmatrix} 1 & 0\\0 & b\end{bmatrix} = \begin{bmatrix} 1 & 0\\0 & b\end{bmatrix} V$$

As $W$ is unimodular, it has integer inverse, so all matrices above are integer. The matrix on the left has second column divisible by $b$, while the matrix on the right has second row divisible by $b$. As they equate, both sides must equal a matrix of the form $\begin{bmatrix} x & by\\bz & bw\end{bmatrix}$ where $w,x,y$ and $z$ are integers. But then $V$ equals $\begin{bmatrix} x & by\\z & w\end{bmatrix}$, so $V$ unimodular implies $xw - byz = \pm 1$.

Similarly, we have $W^{-1} M W = \begin{bmatrix} x & y\\bz & w\end{bmatrix}$, since the diagonal matrix on the right simply scales the second column by $b$. So we started with a general subgroup $H$ invariant under $M$ - generated by the columns of a matrix $A$ - and showed that this $H$ has a basis of the form $W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}$ where the index $b$ divides the bottom left entry of $W^{-1} M W$, as claimed.

For the other direction, let $b$ be a positive factor of the bottom left entry of $W^{-1} M W$, where $W$ is any unimodular matrix. This means

$$W^{-1} M W = \begin{bmatrix} x & y\\bz & w\end{bmatrix}$$

for some integers $w, x, y$ and $z$, where $xw - bzy =\pm 1$, as the matrix is unimodular. But then
$$ M W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix} = W \begin{bmatrix} x & y\\bz & w\end{bmatrix}\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix} $$
and
$$M \left(W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}\right)= W\begin{bmatrix} x & by\\bz & bw\end{bmatrix}=W \begin{bmatrix} 1 & 0\\0 & b\end{bmatrix} \begin{bmatrix} x & by\\z & w\end{bmatrix}= \left(W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}\right) V $$

where the matrix $V$ is unimodular, as we have already noted $xw - bzy =\pm 1$. But this means the subgroup $H$ of $\mathbb{Z}^2$ generated by the columns of $ W\begin{bmatrix} 1 & 0\\0 & b\end{bmatrix}$ is invariant under $M$, as claimed.