Evaluating $\int_0^\infty (\operatorname{E}_n(x)e^x-\frac1{1+x})dx$

I want to evaluate $$I_n=\int_0^\infty \left(\operatorname{E}_n(x)e^x-\frac1{1+x}\right)dx=-\psi(n)$$where $$\operatorname{E}_n$$ denotes exponential integral and $$\psi$$ denotes polygamma function.

Attempt
$$I_n=\int_0^\infty\left(\int_1^\infty\frac{e^{xt}e^x}{t^n}dt-\frac1{1+x}\right)dx\\ =\int_0^\infty\left(\int_0^\infty\frac{e^{xt}}{(t+1)^n}dt-\frac1{1+x}\right)dx\\ =\int_0^\infty\left(\int_0^\infty\frac{e^{xt}}{(t+1)^n}-e^{-(x+1)t}\right)dtdx$$ I tried to change the position of two integral signs but failed because of divergency.

Since $$\operatorname{E}_n (x) \equiv \int \limits_1^\infty \frac{\mathrm{e}^{-x t}}{t^n} \, \mathrm{d} t \sim \frac{\mathrm{e}^{-x}}{x}\left[1 - \frac{n}{x} +\mathcal{O}\left(\frac{1}{x^2}\right)\right]$$ holds as $$x \to \infty$$ , $$I_n$$ exists for $$n \in \mathbb{N}$$ . We have $$\operatorname{E}_n' = - \operatorname{E}_{n-1}$$ for $$n \in \mathbb{N}$$ by definition, so we can integrate by parts to obtain $$I_n = \operatorname{E}_n(x) \mathrm{e}^x ~\Bigg \vert_{x=0}^{x=\infty} + \int \limits_0^\infty \left[\operatorname{E}_{n-1}(x) \mathrm{e}^x - \frac{1}{1+x}\right] \, \mathrm{d} x = - \frac{1}{n-1} + I_{n-1}$$ for $$n > 1$$ . This implies $$I_n = I_1 - H_{n-1}$$ for $$n \in \mathbb{N}$$ . In order to evaluate $$I_1$$ we integrate by parts again, now using the antiderivative $$x \mapsto \mathrm{e}^x - 1$$ of $$x \mapsto \mathrm{e}^x$$ , however. We find $$I_1 = \operatorname{E}_1(x) (\mathrm{e}^x-1) ~\Bigg \vert_{x=0}^{x=\infty} + \int \limits_0^\infty \left[\frac{\mathrm{e}^{-x}}{x} (\mathrm{e}^x - 1) - \frac{1}{1+x}\right] \, \mathrm{d} x = \int \limits_0^\infty \left[\frac{1}{1+x} - \mathrm{e}^{-x}\right] \, \frac{\mathrm{d} x}{x} = \gamma$$ from a well-known integral representation of the Euler-Mascheroni constant, whence the desired result follows: $$I_n = \gamma - H_{n-1} = - \psi(n) \, , \, n \in \mathbb{N} \, .$$