Flipping the gamma function and looking at Newton interpolation provides another angle on the question:
With $\bigtriangledown^{s-1}_{n}c_n=\sum_{n=0}^{\infty}(-1)^n \binom{s-1}{n}c_n$ and $Real(s)>0$, Newton interpolation provides
$\frac{x^{s-1}}{(s-1)!}=\bigtriangledown^{s-1}_{n} \bigtriangledown^{n}_{j} \frac{x^j}{j!}$.
Consistently, interpolating $a_j=\int_0^\infty e^{-x} \frac{x^{j}}{j!}dx=1$ leads to $\int_0^\infty e^{-x} \frac{x^{s-1}}{(s-1)!}dx=1$.