It is well-known applying the monotone convergence theorem that for all $x>1$ we have the following functional equation : $\zeta(x)\Gamma(x)= \displaystyle \int \limits_{[0,+\infty]} \dfrac{t^{x-1}}{\exp(t)-1}\mathrm{d}\lambda(t)$ or $\zeta(x)=\dfrac{1}{\Gamma(x)}\displaystyle \int \limits_{[0,+\infty]} \dfrac{t^{x-1}}{\exp(t)-1}\mathrm{d}\lambda(t) \ $ (where $\lambda$ is the Lebesgue's measure).
Recall that : $\Gamma(x) = \displaystyle \int \limits_{[0,+\infty]}\exp(-t)t^{x-1}\mathrm{d}\lambda(t)$ and $\zeta(x)= \displaystyle \sum\limits_{n=1}^{+\infty} \dfrac{1}{n^x}$ and are well-defined for all $x>1$.
Is it possible to find a link between this functional equation and Bernoulli's numbers ?
I was thinking about the following result (which is not trivial) : $\dfrac{t}{\exp(t)-1} = \displaystyle \sum\limits_{n=0}^{+\infty}B_n \dfrac{t^n}{n!} \ $ where the sequence $(B_n)_{n\ge 0}$ represents the Bernoulli's numbers.
Then we will obtain for all $x > 1$ : $\zeta(x)\Gamma(x)= \displaystyle \int \limits_{[0,1]}t^{x-2}\left( \displaystyle \sum\limits_{n=0}^{+\infty}B_n \dfrac{t^n}{n!} \right) \mathrm{d}\lambda(t) + \displaystyle \int \limits_{[1,+\infty]} \dfrac{t^{x-1}}{\exp(t)-1} \mathrm{d}\lambda(t) \\ = \displaystyle \int \limits_{[0,1]}\left( \displaystyle \sum\limits_{n=0}^{+\infty}B_n \dfrac{t^{n+x-2}}{n!} \right) \mathrm{d}\lambda(t) + \displaystyle \int \limits_{[1,+\infty]} \dfrac{t^{x-1}}{\exp(t)-1} \mathrm{d}\lambda(t) \\ = \displaystyle\sum \limits_{n=0}^{+\infty} \dfrac{B_n}{n!(n+x-1)} + \displaystyle \int \limits_{[1,+\infty]} \dfrac{t^{x-1}}{\exp(t)-1} \mathrm{d}\lambda(t)$
The last equality is obtained using the fact that the series converges uniformly on every compact set included in $[-2\pi, 2\pi]$.
I doubt that we could go further ?
Thnaks in advance !
My own opinion, but still worth posting as an answer, not a comment. The identities $$\zeta(-n)=(-1)^n\frac{B_{n+1}}{n+1},\qquad\zeta(2n)=(-1)^{n+1}\frac{(2\pi)^{2n}}{2(2n)!}B_{2n}\qquad(n\in\mathbb{Z}_{\geqslant0})$$ look like a complete story (along the lines of the question).
With Riemann's functional equation at hand, each one is a consequence of the other; so it's the connection to the (definition of) Bernoulli numbers that deserves a discussion (although I'm sure it's done in a number of other questions here on MSE).
Taking the generating function $z/(e^z-1)=\sum_{n=0}^\infty B_n z^n/n!$ as the definition, the second of the identities is obtained from $$\sum_{n=1}^\infty\zeta(2n)z^{2n}=\sum_{k,n=1}^\infty\left(\frac zk\right)^{2n}=\sum_{k=1}^\infty\frac{z^2}{k^2-z^2}=\frac12(1-\pi z\cot\pi z)$$ using a known consequence of the infinite product for the sine function, and then the complex-exponential form of $\cot\pi z$, making a direct link to the definition of Bernoulli numbers as stated.
An alternative way to get the first of the identities above (i.e., not via the functional equation) is to play Ramanujan's master theorem back, starting with (your well-known representation) $$\Gamma(s)\zeta(s)=\int_0^\infty\frac{x^{s-1}\,dx}{e^x-1}.$$