I also always wondered about the appearance of the quite arbitrary looking function $\dfrac{x}{\exp(x)-1}$. Here is one way of not letting it come out of the blue but getting it by a line of argument. This is an answer to the p.s. of the questioner, so only a partial answer to the initial question, since it will not contribute to the question why the Bernoulli numbers appear everywhere (which, as far as I can see, has not been satisfactorily answered in this thread, anyway), but will refer to the aspect that they are natural things to consider.
As, by the way, are the Bernoulli polynomials, which on our way of developing things also make a natural appearance.
The point of departure is the determination of the power sum polynomial $S_p(x)$, $p = 0,1,2,\dots$. This is the unique polynomial with
$$
\forall{k \in \mathbb{N}:\quad S_p(k) = 0^p+1^p+2^p+\cdots+(k-1)^p},
$$
where $0^p := [p=0]$. It cann alternatively be characterized as the unique polynomial $S_p(x)$ satisfying
$$\tag{1}
\Delta S_p(x) = x^p \quad,\quad S_p(0) = 0.
$$
where $\Delta$ is the forward difference operator acting on polynomials via
$$
\Delta f(x) := f(x+1) - f(x).
$$
Differentiating (1) gives, since $\Delta$ commutes with differentiation,
$$
\Delta S'_p(x) = px^{p-1} = p\Delta S_{p-1}(x)
$$
or
$$
\Delta(S'_p(x) - pS_{p-1}(x)) = 0
$$
or
$$\tag{1a}
S'_p(x) - pS_{p-1}(x) = B_p
$$
for some constant $B_p$, (dependent on $p$). To steer this ship into familiar waters, define
polynomials
$$\tag{2}
B_p(x) := S'_p(x) \quad,\quad p=0,1,2,\dots
$$
of degree $p$ as to have, on differentiating (2):
$$\tag{3}
B'_p(x) - pB_{p-1}(x) = 0,
$$
or
$$\tag{4}
\fbox{$B'_p(x) = pB_{p-1}(x)\quad,\quad p=0,1,2,\dots$}
$$
so that we have done away succesfully with the unwanted (for the time being) constants $B_p$.
The condition $S_p(0)=0$ can be transformed, since $S_p(1) = [p=0]$ by definition of $S_p(x)$, into the equivalent condition
$$\tag{5}
S_p(1) - S_p(0) = \int_0^1 B_p(t) dt = [p=0].
$$
Thus we have reformulated the problem of finding $S_p(x)$ of degree $p+1$ as finding $B_p(x)$ of degree $p$ such that (4) and (5) hold; then we recover $S_p(x)$ as
$$\tag{5a}
S_p(x) = \int_0^x B_p(t) dt \overset{(3)}{=} \dfrac{1}{p+1} \int_0^x B'_{p+1}(t) dt
= \dfrac{1}{p+1}(B_{p+1}(x)-B_{p+1}(0)),
$$
and so we have to take care of the $B_p(x)$, in particular to show that they are well-defined by the conditions (4) and (5).
As said before, the intention behind the introduction of the $B_p(x)$ was to reach familiar waters, and, indeed, we have done so, because (4) shows that the sequence of the $B_p(x)$ is an Appell sequence, so named after Paul Émile Appell (1855-1930):
Definition A sequence of polynomial $A_p(x), p=0,1,2,\dots$ is called an Appell sequence if
(i) $A_0(x)$ is a nonzero constant;
(ii) $\forall{p=1,2,\dots}\quad\quad A'_p(x) = pA_{p-1}(x)$.
Since in an Appell sequence $A_p(x)$ is, by integration, determined by $A_{p-1}(x)$ up to a constant, $a_p$, say, via
$$\tag{6}
A_p(x) = p\int_0^x A_{p-1}(t) dt + a_p,
$$
where $a_0 \not= 0$ because of (i), the Appell sequences over $\mathbb{R}$ -- or any field of characteristic 0 -- correspond in a one-to-one fashion to invertible power series
$$
a(x) = \sum_{p=0}^{\infty} a_p\dfrac{t^p}{p!} \in \mathbb{R}[[t]].
$$
If the Appell sequence $(A_p(x))_{p \in \mathbb{N}}$ corresponds to $a(x)$ via (6), $a(x)$ is called its characteristic series.
To identify the Appell sequence $(B_p(x))_{p \in \mathbb{N}}$, our task is to identify its characteristic series $b(x)$, and for this we need more theory. So consider, for a given Appell sequence $(A_p(x))_{p \in \mathbb{N}}$ its Maclaurin series
$$
A_p(x) = \sum_{k=0}^p A^{(k)}_p(0) \dfrac{x^k}{k!}.
$$
Because of (ii), one has
$$
A^{(k)}_p(0) = p(p-1)\cdots(p-k+1)A_{p-k}(0) = p(p-1)\cdots(p-k+1)a_{p-k},
$$
so that
$$\tag{7}
A_p(x) = \sum_{k=0}^p \binom{p}{k} a_{p-k}x^k.
$$
This is a binomial convolution, so it suggests to introduce the exponential generating function (EGF) of any sequence $(f_p(x))_{p \in \mathbb{N}}$ of polynomials as
$$
f(x,t) := \sum_{p=0}^{\infty} f_p(x) \dfrac{t^p}{p!} \in \mathcal{P}[[t]]
$$
with $\mathcal{P} := \mathbb{R}[t]$ the polynomial ring. For example, if $f(x) = x^p$, then
$$
f(x,t) = \sum_{p=0}^{\infty} x^p \dfrac{t^p}{p!} = \exp(xt).
$$
Then, for an Appell sequence $(A_p(x))_{p \in \mathbb{N}}$, with
$$\tag{8}
A(x,t) := \sum_{p=0}^{\infty} A_p(x) \dfrac{t^p}{p!},
$$
one has
$$\tag{9}
A(0,t) = a(t) \quad,\quad A(x,t) = A(0,t)\exp(xt) = a(t)\exp(xt).
$$
Namely, if I put $x:=0$ in (8):
$$
A(0,t) = \sum_{p=0}^{\infty} A_p(0) \dfrac{t^p}{p!} =:a(t),
$$
which is the first claim of (9), while (7) implies the second claim of (9). These relations show explicitely how the Appell sequence $(A_p(x))_{p \in \mathbb{N}}$ and its characteristic series $a(t)$ determine each other: given $A(x,t)$, one has $a(t) = A(0,x)$, and given $a(x)$, one has $A(x,t) = a(t)\exp(xt)$, so that the $A_p(x)$ are just given by the binomial convolution of the sequences $(a_0,a_1,a_2,\dots)$ and $(1,x,x^2,\dots)$.
By (6), we see that at any stage the polynomial $A_p(x)$ of an Appell sequence is determined by its predecessor up to a constant which is nonzero for $p=0$, so this triggers the idea of fixing an Appell sequence uniquely by putting a nondegenerate linear constraint on each member; in fact, this is exactly what happens in (5), since
$$
f(x) \mapsto \int_0^1 f(t) dt
$$
is a linear form on $\mathcal{P}$. To explicate this, let $L:\mathcal{P} \rightarrow \mathbb{R}$ a linear form on $\mathcal{P}$. We denote the value of $L$ on the polynomial $f(x) \in \mathcal{P}$ by $L[f(x)]$. Let $F:=(f_p(x))_{p \in \mathbb{N}}$ be a basic sequence, i.e. having $\deg f_p(x) = p$ for all $p$. Since a linear form is uniquely defined by the sequence of values it takes on a basic sequence, the linear form $L$ is uniquely defined by, and uniquely defines, the power series
$$
c_{L,F}(t) := \sum_{p=0}^{\infty} L[f_p(x)] \dfrac{t^p}{p!}
$$
In particular, if $M:=(x^p)_{p \in \mathbb{N}}$ is the standard monomial basis, this is the case for
$$
c_L(t) := c_{L,M}(t) = \sum_{p=0}^{\infty} L[x^p] \dfrac{t^p}{p!}
$$
which we call the canonical series of $L$, and we call $L$ nondegenerate if $L[1] \not= 0$ (this terminology is not standard). Then
Theorem Let $L \in \mathcal{P}^*$ be a nondegenerate linear form on $\mathcal{P}$.
Then there is a unique Appell sequence $A:=(A_p(x))_{p \in \mathbb{N}}$ with
$$\tag{10}
\forall p \in \mathbb{N}:\quad L[A_p(x)] = [p=0],
$$
namely the unique Appell sequence with characteristic series
$$\tag{11}
a(t) = (c_{L}(t))^{-1},
$$
i.e. with EGF
$$\tag{12}
A(x,t) = (c_{L}(t))^{-1} \exp(xt).
$$
Proof $\,\,$ By (7) we have, for any Appell sequence $A:=(A_p(x))_{p \in \mathbb{N}}$:
$$\
L[A_p(x)] = \sum_{k=0}^p \binom{p}{k} a_{p-k}L[x^k]
$$
which just says
$$\
c_{L,A} = a(t) c_L(t).
$$
But the conditions (10) are just equivalent to $c_{L,A} = 1$, and so we are done. QED
So, to determine the Appell sequence $B = (B_p(x))_{p \in \mathbb{N}}$, all we have to do is to determine its characteristic series $b(t)$ and so, by (11) and (12), the canonical series $c_L(t)$ for $L$ given by (5):
$$
c_L(t) = \sum_{p=0}^{\infty} \left(\int_0^1 x^p dx \right)\dfrac{t^p}{p!}
= \sum_{p=0}^{\infty} \dfrac{1}{p+1} \dfrac{t^p}{p!}
= \sum_{p=0}^{\infty} \dfrac{t^p}{(p+1)!}
= \dfrac{\exp(t)-1}{t},
$$
which gives us the EGF of its constant coefficients $B_p:=B_p(0)$ as
$$\tag{13}
b(t) = \left(\sum_{p=0}^{\infty}
\dfrac{t^p}{(p+1)!}\right)^{-1} =
\dfrac{t}{\exp(t)-1},
$$
and the EFG of the $B_p(x)$ as
$$\tag{14}
B(x,t) = b(t)\exp(xt) = \dfrac{t\exp(xt)}{\exp(x)-1}.
$$
hence so we have recovered the formerly unwanted constants $B_p$ of (1a) and call them Bernoulli numbers; they are given, but not defined, by their EGF(13). And so have we the polynomials $B_p(x)$ and call them Bernoulli polynomials; they are given, but not defined, by their EGF(14). And, in passing, we have proved, via (5a), Jakob Bernoulli's celebrated formula
$$\
S_p(x) = \dfrac{1}{p+1} \sum_{k=1}^{p+1} \binom{p+1}{k}B_{n-k+1} x^k,
$$
so to speak the cradle of the Bernoulli numbers.
Final Remarks It is an annoying fact that it has become common practice to define the Bernoulli numbers by their EGF (13) (this EGF goes back to Euler, by the way), since this ad hoc definition comes out of the blue and has motivating value approaching zero exponentially. The same holds for the common definition of the Bernoulli polynomials via (14),
which comes staggering along with the same motivational value. The computation above makes them appear, at least, as a result of a comprehensible line of argument. Unfortunately, however, this line of argument appears to add nothing to the clearance of the question why the Bernoulli numbers are so ubiquitous.
