Estimation of the number of solutions by Counting
Consider a binomial random variable $X \sim \mathrm{bin}(n,p)$, and let $Y = (X - pn)^2$. A straightforward but long calculation shows that
$$
\frac{\mathbb{E}[Y]^2}{\mathbb{E}[Y^2]} =
\frac{np(1-p)}{p^3 + (1-p)^3 + 3(n-1)p(1-p)} = \frac{1}{3} + o(1).
$$
(This is to be expected, since a binomial random variable has roughly Gaussian distribution, and for Gaussian distributions the ratio is exactly $1/3$.)

Now let $Z = (S - M)^2$. It is still the case that
$$
\frac{\mathbb{E}[Z]^2}{\mathbb{E}[Z^2]} = \frac{1}{3} + o(1),
$$
since $Z$ is just a scalar multiple of $Y$ above, for an appropriate choice of parameters.

In view of applying the Paley–Zygmund inequality, let us estimate
$$
\theta = \frac{M}{\mathbb{E}[Z]} = \frac{M}{M(N-M)/k} = \frac{k}{N-M}.
$$
It follows that
$$
\Pr[|S-M| &gt; \sqrt{M}] = \Pr[Z &gt; M] \geq (1-\theta)^2 \left(\frac{1}{3} + o(1)\right).
$$
If we want the right-hand side to be bounded by some constant (in your case, $1/4$) then we need $\theta$ to be bounded from below by some constant $c$ (in your case, roughly $1-\sqrt{3/4}$). In other words, we need
$$
k = \theta (N-M) = \Omega(N-M).
$$
As long as $M$ is not very close to $N$, this is $\Omega(N)$.

(When $M$ is very close to $N$, less samples are needed. For example, only one sample is needed when $M = N$.)