Let $\Gamma(z)$ and $\zeta(z)$ denote the gamma and the zeta functions respectively, namely

$$\begin{aligned} &\Gamma(z)=\int_{0}^{\infty} x^{z-1} e^{-x} d x, \quad \operatorname{Re} z>0 \\ &\zeta(z)=\sum_{m=1}^{\infty} \frac{1}{m^{z}}, \quad \operatorname{Re} z>1 \end{aligned}$$

By employing a series expansion of $\left(1-e^{-x}\right)^{-2}$, prove the following identity

$$\int_{0}^{\infty} \frac{x^{z}}{\left(e^{x}-1\right)^{2}} d x=\Gamma(z+1)[\zeta(z)-\zeta(z+1)], \quad \operatorname{Re} z>1$$