By Taylor expanding
$$\frac{x}{e^{x}-1} = \sum_{n=0}^\infty \frac{B_n}{n!}x^n$$
in the Zeta function
$$\zeta(s) = \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} (\frac{x}{e^{x}-1})dx$$
we find
\begin{align} \zeta(s) &= \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} (\frac{x}{e^{x}-1})dx = \frac{1}{\Gamma(s)} \int_0^\infty x^{s-2} \sum_{n=0}^\infty \frac{B_n}{n!}x^n dx \\ &= \frac{1}{\Gamma(s)} \sum_{n=0}^\infty \frac{B_n}{n!} \int_0^\infty x^{s + n - 2} dx \end{align}
what exactly do we do next to arrive at the formulas for $\zeta(-n)$ and $\zeta(2n)$?
By Taylor expanding, for $|x| < 2 \pi$ :
$$\frac{x}{e^{x}-1} = \sum_{n=0}^\infty \frac{B_n}{n!}x^n$$
in the Zeta function, for $Re(s) > 1$ :
$$\zeta(s)\Gamma(s) = \int_0^\infty \frac{x^{s-1}}{e^{x}-1}dx = \int_0^a x^{s-2} \frac{x}{e^{x}-1}dx + \int_a^\infty \frac{x^{s-1}}{e^{x}-1}dx$$
we find whenever $0 < a < 2 \pi$ and $Re(s)> 1$, inverting $\sum$ and $\int$ by absolute/monotone/dominated convergence :
\begin{align} \zeta(s)\Gamma(s) &= \int_a^\infty \frac{x^{s-1}}{e^{x}-1}dx + \int_0^a x^{s-2} \sum_{n=0}^\infty \frac{B_n}{n!}x^n dx \\ &= \int_a^\infty \frac{x^{s-1}}{e^{x}-1}dx + \sum_{n=0}^\infty\frac{B_n}{n!} \int_0^a x^{s + n - 2} dx\\ &= \int_a^\infty \frac{x^{s-1}}{e^{x}-1}dx + \sum_{n=0}^\infty\frac{B_n}{n!} \frac{a^{s+n-1}}{s+n-1}\end{align}
note how by analytic continuation it stays valid for every $s \in \mathbb{C}$ except at the poles,
and how it tells us the residue of $\zeta(s) \Gamma(s)$ at its poles, and hence the value of $\zeta(-k)$ for $k \in \mathbb{N}$