To my surprise, I was able to evaluate the following expression in Mathematica:
$$E\left[e^X \left(1-(1-e^{-X}\right)^n) \right] = \frac{y}{y-1} \left(1-\frac{1}{\binom{n+y-1}{y-1}}\right)\quad X\sim\text{Exp}(y)$$ with the right hand side being equal to $\text{HarmonicNumber}(n)$ in the particular case $y=1$, and equal to $n$ when $y=0$.
If we define $\binom{n}{k} = \frac{\Gamma(n+1)}{\Gamma(n-k+1)\Gamma(k+1)}$ the results seems to hold for all $x,y\in\mathbb R$, though I mostly care about $n$ and $y$ as positive integers.
I have no idea how to prove this by hand. I tried a series expansion of the exponentials without realizing much. I also tried rewriting in terms of the uniform distribution, since $e^{-X}$ has the same distribution as $U^{1/y}$ for $U\sim\text{Uniform}(0,1)$.
Are there any tricks or properties I'm missing?