For k=0,1,..., let ak be the rational number defined by the serie expansion
1ex−1=1x+∑∞k=0(−1)kakΓ(k+1)xk
The computation of the question simply shows that if there would be a "good value" to attach to ∑+∞n=1nk, it should be ak. The real question is why this "good value" is the same as the one obtained via the zeta function: ζ(−k). The answer is that the analytic continuation of the zeta function is defined by essentially the same computation.
The zeta function is defined by the formula ζ(s)=∑n≥1n−s for s∈C, Re(s)>1. To show that this function admits a meromorphic continuation to the all complex plane, a way to do is to use the trick to write
n−s=1Γ(s)∫+∞0e−ntts−1dt.
(Remark not useful for what follows: this trick is used very often in quantum field theory, where t is called in this context a "Schwinger parameter" or the "Schwinger proper time").
Using this and summing the geometric serie, one finds for Re(s)>1.
ζ(s)=1Γ(s)∫+∞0tet−1ts−2dt.
It is possible to do an integration by part to obtain:
−1(s−1)Γ(s)∫+∞0ddt(tet−1)ts−1dt
The point is that this expression makes sense for Re(s)>0 so this formula defines a (and so the) analytic continuation of the zeta function to Re(s)>0. Doing successive integrations by part, one obtains the "full" zeta function and it is immediate from this that ζ(−k)=ak.
This computation defining the analytic continuation of the zeta function is really the same as (precisely the Mellin transform of) the OP's computation. It is one of the standard way to guess what ∑+∞n=1nk "should be".
Remark: essentially the same computation appears in other ways to "regularize" ∑+∞n=1nk, for example the exponential regularization ∑+∞n=1nke−nϵ