$\newcommand{\si}{\sigma}\newcommand{\Si}{\Sigma}$By rescaling, without loss of generality $\lambda=1$. 

The random variable
\begin{equation}
	U_n:=\frac{S-ES}{\binom nk}=\frac1{\binom nk}\sum_{J\in\binom{[n]}k}(Y_J-EY_J)
\end{equation}
is a [U-statistic][3]. Therefore, for each natural $k\ge2$, by Hoeffding's [Theorem 7.1][1], $U_n$ is asymptotically normal (as $n\to\infty$) with (asymptotic) mean $0$ and asymptotic variance $k^2\si_1^2/n$, where 
\begin{equation}
	\si_1^2:=Var\,g(X_1),\quad g(X_1):=E(Y_{[k]}|X_1). 
\end{equation}
In our case, 
\begin{equation}
	g(x)=1(x\ge t)+1(x<t)(1-q^{k-1}),\quad q:=1-e^{-t}, 
\end{equation}
and hence 
\begin{equation}
	\si_1^2=1-q+q \left(1-q^{k-1}\right)^2-\left(1-q^k\right)^2=(1-q)q^{2k-1}. 
\end{equation}

It is easy to see that $ES=\binom nk (1-q^k)$. 

Thus, $S$ is asymptotically normal with (asymptotic) mean $ES=\binom nk (1-q^k)$ and asymptotic variance $\Si^2:=\binom nk ^2 k^2\si_1^2/n$. 

For $n=100$, $k=2$, and $t=1$, we get $S\approx N(ES,\Si^2)$ with 
$ES\approx2972$ and $\Si:=\sqrt{\Si^2}\approx302$, which is in agreement with your picture. 

Explicit bounds on the rate of convergence of $U_n$ to normality (obtained using Stein's method) are available -- see [Chen and Shao, Theorem 3.1][2].

[1]: https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-19/issue-3/A-Class-of-Statistics-with-Asymptotically-Normal-Distribution/10.1214/aoms/1177730196.full

[2]: https://projecteuclid.org/journals/bernoulli/volume-13/issue-2/Normal-approximation-for-nonlinear-statistics-using-a-concentration-inequality-approach/10.3150/07-BEJ5164.full 

[3]: https://en.wikipedia.org/wiki/U-statistic