I have seen Bernoulli numbers many times, and sometimes very surprisingly. They appear in my textbook on complex analysis, in algebraic topology, and of course, number theory. Things like the criteria for regular primes, or their appearance in the Todd class, zeta value at even numbers looks really mysterious for me. (I remember in Milnor's notes about characteristic class there is something on homotopy group that has to do with Bernoulli numbers, too, but I don't recall precisely what that is. I think they also arise in higher K-theory.) 

The list can go on forever. And the wikipedia page of Bernoulli number is already quite long.

My question is, why do they arise everywhere? Are they a natural thing to consider?

==========================================

p.s.----(maybe this should be asked in a separate question) 

Also, I've been wondering why it is defined as the taylor coefficient of the particular function $\frac{x}{e^x-1}$, was this function important? e.g. I could have taken the coefficient of the series that defines the L-genus, namely $\dfrac{\sqrt{z}}{\text{tanh}\sqrt{z}}$, which only amounts to change the Bernoulli numbers by some powers of 2 and some factorial. I guess many similar functions will give you the Bernoulli numbers up to some factor. Why it happen to be the function $\frac{x}{e^x-1}$?

I don't know of a universal theory of all places where Bernoulli numbers arise, but Euler-Maclaurin summation explains many of their more down-to-earth occurrences.

The heuristic explanation (due to Lagrange) is as follows.  The first difference operator defined by $\Delta f(n) = f(n+1)-f(n)$ and summation are inverses, in the same sense in which differentiation and integration are inverses.  This just amounts to a telescoping series: $\sum_{a \le i < b} \Delta f(i) = f(b) - f(a)$.

Now by Taylor's theorem, $f(n+1) = \sum_{k \ge 0} f^{(k)}(n)/k!$ (under suitable hypotheses, of course).  If we let $D$ denote the differentation operator defined by $Df = f'$, and $S$ denote the shift operator defined by $Sf(n) = f(n+1)$, then Taylor's theorem tells us that $S = e^D$.  Thus, because $\Delta = S-1$, we have $\Delta = e^D - 1$.

Now summing amounts to inverting $\Delta$, or equivalently applying $(e^D-1)^{-1}$.  If we expand this in terms of powers of $D$, the coefficients are Bernoulli numbers (divided by factorials).  Because of the singularity at "$D=0$", the initial term involves antidifferentiation $D^{-1}$, i.e., integration.  Thus, we have expanded a sum as an integral plus correction terms involving higher derivatives, with Bernoulli number coefficients.

Specifically,
$$
\sum_{a \le i < b} f(i) = \int_a^b f(x) \, dx + \sum_{k \ge 1} \frac{B_k}{k!} (f^{(k-1)}(b) - f^{(k-1)}(a)).
$$
(Subtracting the values at $b$ and $a$ just amounts to the analogue of turning an indefinite integral into a definite integral.)

This equation isn't literally true in general: the infinite sum usually won't converge and there's a missing error term.  However, it is true when $f$ is a polynomial, and one can bootstrap from this case to the general one using the Peano kernel trick.

So from this perspective, the reason why $t/(e^t-1)$ is a natural generating function to consider is that we sometimes want to invert $e^t-1$ (the factor of $t$ is just to make it holomorphic), and the most important reason I know of to invert it is that we want to invert $\Delta = e^D-1$.
