Coupling showing that $\operatorname{Bin}(n,\frac{1}{n+1})$ stochastically dominates $\operatorname{Bin}(n-1,\frac{1}{n})$ The classical inequality
$$
\left(1-\frac{1}{n}\right)^{n-1} > \frac{1}{e}
$$
has a probabilistic generalization: the binomial distribution $\operatorname{Bin}(n-1,\frac{1}{n})$ is stochastically dominated by the Poisson distribution $\operatorname{Po}(1)$. A simple coupling proof of this can be found in Klenke and Mattner, Stochastic ordering of classical discrete distributions. They also prove the stronger result that $\operatorname{Bin}(n,\frac{1}{n+1})$ stochastically dominates $\operatorname{Bin}(n-1,\frac{1}{n})$, which generalizes the fact that the sequence $(1-\frac1n)^{n-1}$ is decreasing.

Is there a natural coupling that shows that $\operatorname{Bin}(n,\frac{1}{n+1})$ stochastically dominates $\operatorname{Bin}(n-1,\frac{1}{n})$?

What I have in mind is an elementary argument that involves a balls-into-bins scenario, but any other simple coupling would be great.
(Just to be clear: I already know that there exists some coupling, this follows from the stochastic domination.)
 A: I am not sure wether this is a good answer, but I would say monotone coupling would be one possibility. 
Illustration for $n=2$
For example let $n=2$, $\ X\sim \mu=\text{Bin}\bigl(2, \frac13 \bigr)$, $ \ Y\sim \nu=\text{Bin}\bigl(1,\frac12 \bigr)$.
Since $X$ can only take three values a natural candidate for a suitable probability space for $X$ would be $\Omega = \{\omega_0, \omega_1, \omega_2\}$
where $X(\omega_i)=i$.
But it turns out that we have to split $\omega_1$ into two parts to achieve monotone coupling. So we pick $\Omega = \{\omega_0, \omega_{1a}, \omega_{1b}, \omega_2\}$ and define $X$ as $X(\omega_{i})=i$ for $i=0,2$ and $X(\omega_{1a}) = X(\omega_{1b})=1$. This will become clear later.
Then we define a probability measure $\mathbb{P}$ such that $X\sim \mu$ under $\mathbb{P}$, namely 
$$
\mathbb{P}(\omega_0) = \frac49, \quad \mathbb{P}(\omega_{1a}) = \frac1 {18}, \quad 
\mathbb{P}(\omega_{1b}) = \frac 7 {18}, \quad \mathbb{P}(\omega_{2}) = \frac19.
$$
We want to construct now a random variable $Y$ on $\Omega$ such that $Y \sim \nu$ under $\mathbb{P}$ and we do this by monotone coupling. 
The mass of $X$ (or resp $\mu$) is monotonly transferred to $Y$ (or respectively to $\nu$).
In the end $Y$ should only take two values (0 and 1) with probability $\frac12$ each.
We start from the left and put all of $\mu$'s mass from zero (i.e. $\omega_0$) and put it into $\{Y=0 \}$, which means that $Y(\omega_0)=0$. Then there is no more mass left at   $\{X=0\}$ but there is still some room at $\{Y=0 \}$ (since $\mathbb{P}(\omega_0)= \frac 49$ but $\mathbb{P}(Y=0)$ should be $\frac 12$).
So we continue monotonically: we take some mass form $\{X=1\}$ and put it into $\{Y=0\}$ until $\{Y=0\}$ is full (i.e has mass $\frac 12$). Thats why we had to split up $\omega_1$ into $\omega_{1a}$ and $\omega_{1b}$: we can now assign  $\omega_{1a}$ with $\{Y=0 \}$. $\omega_{1a}$ has exactly the right mass to guarantee that $\mathbb{P}(Y=0)= \frac 12$.
The remaining mass at $\{X=1\}$ (i.e $\omega_{1b}$) is transferred to $\{Y=1 \}$, i.e $Y(\omega_{1b})=1$. Then we continue monotonically and transfer the mass from $\{X=2\}$ (which is all the remaining mass) to $\{Y=1 \}$, which means that $Y(\omega_2)=1$.
To sum up
$$
(X,Y)(\omega_0)=(0,0), \quad (X,Y)(\omega_{1a})=(1,0), \quad (X,Y)(\omega_{1b})=(1,1), \quad (X,Y)(\omega_2)=(2,1),
$$
and $X \geq Y$.
General Case
Of course this was only for $n=2$, but the general case is similar. You use the same technique of monotone coupling. Therefore you take all the mass from $\{X=0\}$ and assign it to $\{Y=0\}$. Then you split the mass at $\{X=1\}$, assign one part of it to  $\{Y=0\}$ and the other part to $\{Y=1\}$, and so on.
