Looking for a difference that makes a difference.
Flipping the gamma function and looking at Newton interpolation provides another angle on the question:
Accepting the general binomial coefficient as a natural extension of the integral coefficient through the Taylor series, or binomial theorem, for $(1+x)^{s-1}$ and with the finite differences $\bigtriangledown^{s-1}_{n}c_n=\sum_{n=0}^{\infty}(-1)^n \binom{s-1}{n}c_n$, Newton interpolation gives $$\bigtriangledown^{s-1}_{n} \bigtriangledown^{n}_{j} \frac{x^j}{j!}=\frac{x^{s-1}}{(s-1)!}$$
for $Real(s)>0,$ so the values of the factorial at the nonnegative integers determine uniquely the generalized factorial, or gamma function, in the right-half of the complex plane through Newton interpolation.
Then with the sequence $a_j=1=\int_0^\infty e^{-x} \frac{x^{j}}{j!}dx$,
$\bigtriangledown^{s-1}_{n} \bigtriangledown^{n}_{j}a_j=\bigtriangledown^{s-1}_{n} \bigtriangledown^{n}_{j}1=1=\bigtriangledown^{s-1}_{n} \bigtriangledown^{n}_{j}\int_0^\infty e^{-x} \frac{x^{j}}{j!}dx=\int_0^\infty e^{-x} \frac{x^{s-1}}{(s-1)!}dx$.
This last integral allows the interpolation of the gamma function to be analytically continued to the left half-plane as in MSE-Q132727, so the factorial can be uniquely extended from it's values at the non-negative integers to the entire complex domain as a meromorphic function.
This generalizes in a natural way; the integral, a modified Mellin transform, can be regarded as a means to interpolate the coefficients of an exponential generating function (in this particular case, exp(x)) and extrapolate the interpolation to the whole complex plane. (See my examples in MO-Q79868 and MSE-Q32692 and try the same with the exponential generating function of the Bernoulli numbers to obtain an interpolation to the Riemann zeta function.)
From these perspectives--its dual roles as a natural interpolation of a sequence itself and in interpolating others and also as an iconic meromorphic function--the gamma function provides the "best" generalization of the factorial from the nonneg integers to the complex plane.
To more sharply connect pbrooks interest in fractional calculus and the gamma function with Quiaochu's in the beta integral, it's better to look at a natural sinc interpolation of the binomial coefficient:
Consider the fractional integro-derivative
$\displaystyle\frac{d^{\beta}}{dx^\beta}\frac{x^{\alpha}}{\alpha!}=FP\frac{1}{2\pi i}\displaystyle\oint_{|z-x|=|x|}\frac{z^{\alpha}}{\alpha!}\frac{\beta!}{(z-x)^{\beta+1}}dz=FP\displaystyle\int_{0}^{x}\frac{z^{\alpha}}{\alpha!}\frac{(x-z)^{-\beta-1}}{(-\beta-1)!} dz$
$=\displaystyle\frac{x^{\alpha-\beta}}{(\alpha-\beta)!}$
where FP denotes a Hadamard-type finite part, $x>0$, and $\alpha$ and $\beta$ are real.
For $\alpha>0$ and $\beta<0$, the finite part is not required for the beta integral, and it can be written as
$\displaystyle\int_{0}^{1}\frac{(1-z)^{\alpha}}{\alpha!}\frac{z^{-\beta-1}}{(-\beta-1)!} dz=\sum_{n=0}^{\infty } (-1)^n \binom{\alpha }{n}\frac{1}{n-\beta}\frac{1}{\alpha!}\frac{1}{(-\beta-1)!}$
$=\displaystyle\sum_{n=0}^{\infty }\frac{\beta!}{(\alpha-n)! n!}\frac{\sin (\pi (\beta -n))}{\pi (\beta -n)}=\frac{1}{(\alpha-\beta)!}, $ or
$$\displaystyle\sum_{n=0}^{\infty }\frac{1}{(\alpha-n)! n!}\frac{\sin (\pi (\beta -n))}{\pi (\beta -n)}=\frac{1}{(\alpha-\beta)!\beta!},$$
where use has been made of $\frac{\sin (\pi \beta)}{\pi \beta}=\frac{1}{\beta!(-\beta)!}$, and $\alpha$, of course, can be a positive integer. The final sinc fct. interpolation holds for $Real(\alpha)>-1$ and all complex $\beta$.
Euler's motivation (update July 2014):
R. Hilfer on pg. 18 of "Threefold Introduction to Fractional Derivatives" states, "Derivatives of non-integer (fractional) order motivated Euler to introduce the Gamma function ...." Euler introduced in the same reference given by Hilfer essentially
$\displaystyle\frac{d^{\beta}}{dx^\beta}\frac{x^{\alpha}}{\alpha!}=\frac{x^{\alpha-\beta}}{(\alpha-\beta)!}$.
And, (added Feb. 17,2022), from Whittaker, E., 1928/1929, Oliver Heaviside, The Bulletin of the Calcutta Mathematical Society 20: 199-220:
This [fractional differentiation] is an old subject: Leibniz considered it in 1695 and Euler in 1729: and indeed it was in order to generalize the equations
$\frac{d^n(x^k)}{dx^n}= k(k-1)...(k-n+1)x^{k-n}$ to fractional values of $n$ that Euler invented the Gamma-Function.
Another Useful Interpolation (Edit 1/22/21)
There is a second interpolation of the inverse factorial that nature (and contemporary researchers) seems to find useful involving the Mittag-Leffler function. This again is related to the fractional calculus
From the Cauchy residue theorem, we can represent differentiation via
$$k!\; a(k) = k! \; \oint_{|z|=r} \frac{e^z}{z^{k+1}} \; dz = e^{-1}k! \; \oint_{{|z|=r}} \frac{e^{z+1}}{z^{k+1}} \; dz $$
$$= e^{-1}k! \; \oint_{|z-1|=1} \frac{e^{z}}{(z-1)^{k+1}} \; dz =e^{-1} D^k_{z=1} e^z.$$
Interpolating using a standard fractional integroderivative, of which there are several reps,
$$\lambda! \; a(\lambda) = \; e^{-1} D_{z=1}^{\lambda} \; e^z = e^{-1} \; \sum_{n \ge 0} \frac{z^{n-\lambda}}{(n-\lambda)!}\; |_{z=1}$$
$$ = e^{-z} \; \sum_{n \ge 0} \frac{z^{n-\lambda}}{(n-\lambda)!}\; |_{z=1} = e^{-z} z^{-\lambda}\; E_{1,-\lambda}(z) \; |_{z=1},$$
where $E_{\alpha,\beta}(z)$ is the Mittag-Leffler function (general definition in Wikipedia, MathWorld; some applications), encountered very early on by anyone exploring fractional calculus.
This method of interpolation gives the entire function (over complex $\lambda$)
$$ a(\lambda) = e^{-1} \; E_{1,-\lambda}(1) \frac{1}{\lambda!} = e^{-1} \; \sum_{n \ge 0} \frac{1}{(n-\lambda)!} \; \frac{1}{\lambda!}, $$
which gives $a(k) = \frac{1}{k!} = \frac{1}{\Gamma(k+1)}$ for $k=..., -2,-1,0,1,2, ...$, i.e., the integers.