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?
Denote $$I = E\left[e^X \left(1-(1-e^{-X}\right)^n) \right] =\int_0^{+\infty}e^x(1-(1-e^{-x})^n)ye^{-yx}dx$$ Let's make a change of variables $z = 1-e^{-x}$ then $x = -\ln(1-z)$ and $dx = \frac{1}{1-z}dz$. The integral is from $z = 0$ to $z=1$. We have:
\begin{align} I & = \int_0^{+\infty}y(1-(1-e^{-x})^n)e^{-(y-1)x}dz \\ & = \int_0^{1}y(1-z^n)(1-z)^{y-1}\frac{1}{1-z}dz \\ & = \int_0^{1}y(1-z)^{y-1}(\sum_{k=0}^{n-1} z^k)dz \\ & = \sum_{k=0}^{n-1} \left( \int_0^{1}y(1-z)^{y-1}z^kdz \right) \\ \end{align}
By using Mathematica, we have $$\int_0^{1}y(1-z)^{y-1}z^kdz = \frac{\Gamma(1+k)\Gamma(y+1)}{\Gamma(1+k+y)}$$ Hence, $$I = \sum_{k=0}^{n-1} \left( \frac{\Gamma(1+k)\Gamma(y+1)}{\Gamma(1+k+y)} \right)$$
I resolve the case you asked where $y\in \mathbb{N}^*$, we have
\begin{align} I &= \sum_{k=0}^{n-1} \frac{y!k!}{(k+y)!} \\ &= y!\sum_{k=0}^{n-1} \frac{1}{(k+1)...(k+y)} \\ &= \frac{y!}{y-1}\sum_{k=0}^{n-1} \left( \frac{1}{(k+1)...(k+y-1)} - \frac{1}{(k+2)...(k+y)} \right) \\ &= \frac{y!}{y-1} \left( \frac{1}{(y-1)!} - \frac{n!}{(n+y-1)!} \right) \\ &=\frac{y}{y-1}\left(1-\frac{1}{\binom{n+y-1}{y-1}}\right) \end{align}