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$.
So, is the gamma fct. unique in satisfying the Newton interpolation relation above, which itself is defined in terms of the gamma fct.? I.e., can I substitute another definition of the gamma fct. in the binomial coefficient that satisfies the relation? Any substitute would lead to an unnatural interpolation of the constant $a_j=1$.