How was the integral for Zeta Function created

How was the zeta function integrated from

$$\zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^{s}}$$

To

$$\zeta(s) = \frac{1}{\Gamma (s)}\int_{0}^{\infty}\frac{x^{s-1}}{e^{x}-1}dx$$

I've tried googling this and surprisingly can't find much on it, even the wikipedia article for zeta function doesn't explain how this integral is derived or cite a source anywhere. I have no doubt it's true; I am just curious how it was obtained. I know very little about converting infinite sums to integrals.

Instead of using the dominated or monotone convergence theorem, I like to prove it by elementary means.

For $Re(s) > 0$ and $n > 0$ (change of variable $y = nx$) : $$\Gamma(s)n^{-s} = \int_0^\infty x^{s-1} e^{-nx}dx$$

So that for $Re(s) > 0$ (using the geometric series) : $$\Gamma(s) \sum_{n=1}^N n^{-s} = \int_0^\infty x^{s-1} \sum_{n=1}^N e^{-nx}dx = \int_0^\infty x^{s-1}\frac{1-e^{-Nx}}{e^{x}-1} dx$$

For $Re(s) > 1$ it is known that

$$\zeta(s) = \lim_{N \to \infty} \sum_{n=1}^N n^{-s}$$

Finally, we need to prove that (again for $Re(s) > 1$) : $$\lim_{N \to \infty}\int_0^\infty x^{s-1}\frac{1-e^{-Nx}}{e^{x}-1} dx = \int_0^\infty \frac{x^{s-1}}{e^{x}-1} dx$$ which is obvious once we showed that for $x >0$ : $\displaystyle\left|\frac{x}{e^{x}-1}\right| < 1$ whence $\displaystyle\int_0^\infty \frac{x^{s-1}}{e^x-1}e^{-Nx}dx$ converges absolutely and $\to 0$ as $N \to \infty$

Overall, for $Re(s) > 1$ :

$$\Gamma(s) \zeta(s) = \lim_{N \to \infty} \Gamma(s) \sum_{n=1}^N n^{-s} = \lim_{N \to \infty} \int_0^\infty x^{s-1} \frac{1-e^{-Nx}}{e^x-1}dx = \int_0^\infty \frac{x^{s-1}}{e^{x}-1} dx$$