$\newcommand{\la}{\lambda} \newcommand{\Si}{\Sigma} \newcommand{\thh}{\theta} \newcommand{\vpi}{\varphi} \renewcommand{\P}{\operatorname{\mathsf P}} \newcommand{\E}{\operatorname{\mathsf E}} \newcommand{\Var}{\operatorname{\mathsf Var}}$
Let $X_i:=M[i]$ and $n:=|M|$, so that \begin{equation} X=\sum_{1\le i<j\le n}1_{\{X_i=X_j\}}=\binom n2 U,\quad U:=\binom n2^{-1}X=\binom n2^{-1}\sum_{1\le i<j\le n}h(X_i,X_j), \end{equation} $h(x_1,x_2):=1_{\{x_1=x_2\}}$, so that $U$ is a $U$-statistic with kernel $h$ of order $m=2=\;$the number of arguments of $h$.
We have \begin{equation} \E X=\sum_{1\le i<j\le n}\P(X_i=X_j)=\binom n2\frac1N, \end{equation} and \begin{equation} \E X^2=\sum_{\{i,j\}\in\binom{[n]}2}\;\sum_{\{k,\ell\}\in\binom{[n]}2}\P(X_i=X_j,X_k=X_\ell)=\Si_1+\dots+\Si_4, \end{equation} \begin{equation} \Si_1:=\sum_{\{i,j\}\cap\{k,\ell\}=\emptyset}\P(X_i=X_j,X_k=X_\ell)=\frac1{N^2}\,\binom n2\binom{n-2}2, \end{equation} \begin{equation} \Si_2:=\sum_{\{i,j\}\cap\{k,\ell\}=\{i\wedge j\}}\P(X_i=X_j=X_k)=\frac1{N^2}\,\binom n2\binom{n-2}1, \end{equation} \begin{equation} \Si_3:=\sum_{\{i,j\}\cap\{k,\ell\}=\{i\vee j\}}\P(X_i=X_j=X_k)=\frac1{N^2}\,\binom n2\binom{n-2}1, \end{equation} \begin{equation} \Si_4:=\sum_{\{i,j\}=\{k,\ell\}}\P(X_i=X_j)=\frac1{N}\,\binom n2, \end{equation} whence \begin{equation} \Var X=\E X^2-\E^2X=\frac1N\,\binom n2-\frac1{N^2}\binom{n-1}1. \end{equation}
Now one can use the Paley--Zygmund inequality: for $\thh\in[0,1)$, \begin{multline} \P(X>\thh\E X)\ge 1-\frac{\Var X}{\Var X+(1-\thh)^2\E^2X} \\ = 1-\bigg[N\binom n2-\binom{n-1}1\bigg]\bigg/\bigg[N\binom n2-\binom{n-1}1+(1-\thh)^2\binom n2^2\bigg]\to1 \end{multline} as $n\to\infty$.
The $U$-statistic $U=\binom n2^{-1}X$ is degenerate, since $\Var\E(h(X_1,X_2)|X_1)=\Var\frac1N=0$. Using an appropriate limit theorem for such statistics (see e.g. Theorem 4), we see that $\frac2{n-1}(X-\E X)=n(U-\E U)$ converges in distribution (as $n\to\infty$) to the random variable (r.v.) \begin{equation} \sum_1^N\la_i(Z_j^2-1)=(2/N-1)(Z_1^2-1)+(2/N)\sum_2^N(Z_j^2-1), \end{equation} where the $Z_j$'s are iid standard normal r.v.'s and the $\la_j$ are the eigenvalues of the centered kernel $h(x_1,x_2)-\E h(X_1,X_2)=1_{\{x_1=x_2\}}-\P(X_1=X_2)=1_{\{x_1=x_2\}}-1/N$, so that $\la_1=2/N-1$ (with a corresponding eigenfunction $\vpi_1(x)\equiv1$) and $\la_2=\dots=\la_N=2/N$ with $n-1$ linearly independent eigenfunctions $\vpi_2,\dots,\vpi_N$ satisfying the condition $\sum_{x=1}^N\vpi_j(x)=0$ for $j=2,\dots,N$.