$\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][1] 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][2]: 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][1]), 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$. 

  [1]: https://www.math.ucla.edu/~tom/Stat200C/Ustat.pdf
  
  [2]: https://en.wikipedia.org/wiki/Paley%E2%80%93Zygmund_inequality