According to Wikipedia (https://en.wikipedia.org/wiki/Bernoulli_number) we have: $$1^m + 2^m + ... n^m = S_m(n) = \frac{1}{m+1} \sum_{k=0}^{m} \binom{m+1}{k} B_k^{+} n^{m+1-k}$$ where $B_k^{+}$ are the (positive) Bernoulli numbers.
From plugging in a few values of natural numbers m it seems that $$\int_{-1}^{0} S_m(n) dn = -\frac{1}{m+1} \sum_{k=0}^{m} \frac{1}{m+2-k} \binom{m+1}{k} B_k^{+} (-1)^{m+2-k} = \zeta(-m)$$ where the last function is the Riemann zeta function. Is there any intuitive reason for this?
