Estimate of Bernoulli numbers from the contour integral One knows (see https://mathworld.wolfram.com/BernoulliNumber.html)
that the Bernoulli number $B_n$ is $B_n=\frac{n!}{2i\pi}\int_{\mathcal C}\frac z{e^z-1}\frac{\mathrm dz}{z^{n+1}}$, where $\mathcal C$ is a closed contour included in $\{z\in\mathbb C\mid |z|<2\pi\}$. With this formula, is it possible to obtain the estimate: $|B_{2n}|\sim4\sqrt{\pi n}\left(\frac n{\pi e}\right)^n$? The proof I know depends on the zeta function, that I do not want to use.

Note that the estimate is supposed to be
$\lvert B_{2m} \rvert \sim 4\sqrt{\pi m}\left(\frac{m}{\pi e}\right)^{2m}$, with exponent $2m$, not $m$. And by Stirling's approximation $n! \sim \sqrt{2 \pi n} \left(\frac{n}{e}\right)^n$, we can rewrite it as $\frac{\lvert B_{2m} \rvert}{(2m)!} \sim 2 \left(\frac{1}{2 \pi}\right)^{2m}$.
The contour integral you gave is equivalent to $\frac{z}{e^z-1} = \sum_{n=0}^\infty \frac{B_n}{n!} z^n $. We'll be applying a typical sort of method in combinatorics and number theory, relating asymptotics of a sequence to the singularities of a generating function. The radius of convergence is $2 \pi$ since the singularities nearest to the origin are at $z = \pm 2 \pi i$. From the standard theorem about the radius of convergence of power series, we can already conclude immediately that
$$
  \frac{1}{2 \pi} = \limsup_n \sqrt[n]{\lvert B_n \rvert / n!}.
$$
This is not as strong as the desired result, unfortunately - we need something better.
To improve the previous estimate, consider the function left over after taking $\frac{z}{e^z-1}$ and subtracting off the poles at $z = \pm 2 \pi i$, which are simple poles with residues $\pm 2 \pi i$. That is, consider $\frac{z}{e^z-1} - \frac{2 \pi i}{z - 2 \pi i} - \frac{-2 \pi i}{z + 2 \pi i}$. This has the power series $\sum_{n=0}^\infty \left( \frac{B_n}{n!} + \frac{1}{(2 \pi i)^n} + \frac{1}{(-2 \pi i)^n} \right) z^n$. And now the singularities nearest to the origin are $z = \pm 4 \pi i$. So
$$
  \frac{1}{4 \pi} = \limsup_n \sqrt[n]{\left\lvert \frac{B_n}{n!} + \frac{1}{(2 \pi i)^n} + \frac{1}{(-2 \pi i)^n} \right\rvert}.
$$
In particular, $\frac{B_{2m}}{(2m)!} + (-1)^m \frac{2}{(2 \pi)^{2m}} = O \left( \left( \frac{1}{4 \pi} + \varepsilon \right)^{2m} \right)$ for every $\varepsilon > 0$. This is a stronger form of the desired result $\frac{\lvert B_{2m} \rvert}{(2m)!} \sim 2 \left(\frac{1}{2 \pi}\right)^{2m}$.
