Here's a proof of the positivity of
$$
c_n(\alpha) := \sum_{r=0}^n (-1)^r {n\choose r}^\alpha
$$
for all even $n$ and real $\alpha < 1$. It follows
(via M.Wildon's clever $F(x) F(-x)$ trick at mo.84958) that
$\sum_{n=0}^\infty \phantom. x^n / n!^{\alpha} > 0$ for all $x \in\bf R$.
[EDIT fedja has meanwhile provided a very nice direct proof of
the positivity of $\sum_{n=0}^\infty \phantom. x^n / n!^{\alpha}$.]
The key is to write $c_n(\alpha)$ as a finite difference
$$
\sum_{r=0}^n \phantom. (-1)^r {n\choose r} \cdot {n\choose r}^{\alpha - 1}
$$
and show that the Gamma interpolation
$$
\bigl(\Gamma(r+1)\Gamma(n-r+1) / n!\bigr)^{1-\alpha}
= n!^{\alpha-1} \exp\bigl((1-\alpha) (\log\Gamma(r+1) + \Gamma(n-r+1)\bigr)
$$
of ${n\choose r}^{\alpha - 1}$ has a positive $n$-th derivative
for all $r \in [0,n]$.
This in turn follows from the fact that the expansion of
$\log\Gamma(r+1) + \log\Gamma(n-r+1)$ in a Taylor series about $r = n/2$
has positive $(r - (n/2))^k$ coefficient for each $k=2,4,6,\ldots$.
[The coefficient vanishes for odd $k$ because
$\log\Gamma(r+1) + \log\Gamma(n-r+1)$ is an even function of $r-(n/2)$.]
Indeed the well-known formula
$$
\log \Gamma(x) = -\gamma x - \log x + \sum_{j=1}^\infty
\left[ \frac{x}{j} - \log \left( 1 + \frac{x}{j} \right) \right]
$$
shows that the $k$-th derivative of $\log\Gamma(x)$ is positive
for all $x>0$ and $k=2,4,6,\ldots$, because this is true for
$-\gamma x - \log x$ and for each term in the sum; explicitly
the derivative is $k! \phantom. \sum_{j=0}^\infty (x+j)^{-k}$ which is
positive termwise. Therefore in the Taylor expansion
$$
\log \Gamma(r+1) = \log(n/2)! + \sum_{k=1}^\infty \phantom. g_k (r-(n/2))^k
$$
each of $g_2,g_4,g_6,\ldots$ is even.
Since $\log\Gamma(r+1) + \log\Gamma(n-r+1)$ is
$$
2\log(n/2)!
+ 2 \Bigl( g_2 (r-(n/2))^2 + g_4 (r-(n/2))^4 + g_6 (r-(n/2))^6 + \cdots\Bigr),
$$
the claim follows. [EDIT David Speyer notes that the convergence
of the Taylor series on $|r-(n/2)| \leq n/2$ requires justification,
and that the justification is easy because the $\Gamma(z)$ has no zeros
and poles only at $0,-1,-2,\ldots$ so the radius of convergence is $(n/2)+1$.]
Multiplying by $1 - \alpha$ and substituting into the exponential series,
we deduce that $(\Gamma(r+1) \Gamma(n-r+1))^{1-\alpha}$, too,
is a positive combination of even powers of $r-(n/2)$.
Now if a function $g$ has positive $n$-th derivative, then its
first finite difference
$$
g(x+1) - g(x) = \int_x^{x+1} g'(y) dy
$$
has positive $(n-1)$-st derivative; repeating this argument $n$ times,
we find that the $n$-th finite difference is positive, and we're done.
I played around with your sum in Maple and got
$$ \frac{2n}{n+1}{2n-1 \choose n-1}{n^2 \choose p-n} 3F_{2}([n,n-p,2n+1],[n+2,n^2+n+1-p],1) $$
I make no guarantees that this is correct (especially as the original answer contained a $\binom{n^2}{-1}$ in it).
Best Answer
I have obtained a formula for the generating function of your sequence.
Let $S_n$ be defined as in the quesion. We extend the definition to $n = 0$ by demanding $0^0 = 0$, hence $S_0 = 0$.
Consider $S(t) = \sum_{n\geq 0} S_n t^n$. I will work out a formula for $S(t)$.
\begin{eqnarray*} S(t) &=& \sum_{n=0}^\infty \sum_{j = 0}^n\frac{(-e)^{-j}}{j!}(n - j)^j t^n\\ &=& \sum_{j = 0}^\infty \frac{(-e)^{-j}}{j!}t^j\sum_{n = 0}^\infty n^jt^n\\ &=& \sum_{j = 0}^\infty \frac{(-e^{-1}t)^j}{j!}\frac{tA_j(t)}{(1 - t)^{j + 1}}\\ &=& \frac{t}{1 - t}\sum_{j = 0}^\infty A_j(t)\frac{(\frac{e^{-1}t}{t - 1})^j}{j!}\\ &=& \frac{t}{1 - t}\frac{t - 1}{t - e^{e^{-1}t}}\\ &=& \frac{t}{e^{e^{-1}t} - t}. \end{eqnarray*}
Explanation for the calculation:
The key step is $\sum_{n\geq 0} n^jt^n = \frac{tA_j(t)}{(1 - t)^{j + 1}}$, where $A_j(t)$ is the Eulerian polynomial. We then use the generating function $\sum_{j \geq 0}A_j(t)\frac{x^j}{j!} = \frac{t - 1}{t - e^{(t - 1)x}}$.
It is then reasonable to make the change of variable $T = e^{-1}t$, which leads to $S_n = e^{-n}S'_n$, where the sequence $(S'_n)_n$ has generating function $\frac{eT}{e^T - eT}$. This at least gives a numerically better way to calculate the numbers $S_n$.
And the question becomes to prove that $S'_n = 2n + \frac{2}{3} + o(1)$, where $S'_n$ is the $n$-th coefficient of the Taylor expansion of the function $\frac{T}{e^{T - 1} - T}$ at $T = 0$.
I'm however not able to proceed further...
One may try to subtract $\sum_{n \geq 0}(2n + \frac{2}{3})T^n$, which is a rational fraction, and try to show that the difference has Taylor coefficients tending to $0$. But this doesn't seem to help much.
Also the $n$-th Taylor coefficient could be calculated via a residue at $0$, but again doesn't seem to help much.
Another observation is that $T = 1$ is a pole of the function $\frac{T}{e^{T - 1} - T}$.
Experimentally, it seems that the parameter $e$ is quite important. For this parameter, we indeed have $S'_n = 2n + \frac{2}{3} + o(1)$, as Henri Cohen stated in the comment. In fact, the $o(1)$ is also exponentially decreasing.
If $e$ is changed to anything larger, then the sequence $(S'_n)_n$ increases much faster; if it is changed to anything smaller, then negative terms appear.