A Bernoulli number identity and evaluating $\zeta$ at even integers

As far as I can tell, you have already proven your claim; you started with a function with the appropriate value for the coefficient of interest and then found an expression for that coefficient. Consequently that expression must be a Bernoulli number by construction. If for some reason you are after an alternative proof that the coefficient of $z^{2n}$ in the function:

$g(z)=2(-1)^n(2n)!2^{-2n}f(z)$,

where:

$f \left( z \right) =\prod _{k=0}^{n-1}\mathrm{sinc} \left( z \exp \left( -\dfrac{k \pi i}{n} \right)\right) $,

$\mathrm{sinc}(x)=\dfrac{\sin(x)}{x},$

is equal to the Bernoulli number $B(2n)$ then here's one...

First take the logarithm of $f(z)$:

$\ln(f \left( z \right)) =\sum _{k=0}^{n-1}\ln\left(\mathrm{sinc} \left( z \exp \left(- \dfrac{k \pi i}{n} \right)\right)\right)$.

Then use the product expansion of the sinc function:

$\mathrm{sinc} \left( z \exp \left( -\dfrac{k \pi i}{n} \right)\right)=\prod _{q=1}^{\infty}\left(1-\left(\dfrac{z}{{\pi}q}\right)^2\exp \left(- \dfrac{2k \pi i}{n} \right)\right)$

to obtain:

$\ln(f \left( z \right)) =\sum _{k=0}^{n-1}\sum _{q=1}^{\infty}\ln\left(1-\left(\dfrac{z}{{\pi}q}\right)^2\exp \left(- \dfrac{2k \pi i}{n} \right)\right)$.

then the series expansion of the logarithm:

$\ln\left(1-\left(\dfrac{z}{{\pi}q}\right)^2\exp \left(- \dfrac{2k \pi i}{n} \right)\right)=-\sum _{r=1}^{\infty}\left(\dfrac{z}{{\pi}q}\right)^{2r}r^{-1} \exp \left(- \dfrac{2kr \pi i}{n} \right)$

to get:

$\ln(f \left( z \right)) =-\sum _{k=0}^{n-1}\sum _{q=1}^{\infty}\sum _{r=1}^{\infty}\left(\dfrac{z}{{\pi}q}\right)^{2r}r^{-1} \exp \left(- \dfrac{2kr \pi i}{n} \right)$.

$\ln(f \left( z \right)) =-\sum _{q=1}^{\infty}\sum _{r=1}^{\infty}\left(\dfrac{z}{{\pi}q}\right)^{2r}r^{-1} \sum _{k=0}^{n-1}\exp \left(- \dfrac{2kr \pi i}{n} \right)$.

Note the Kronecker delta popping up that will only pick out $r$ when it's a multiple (denoted $l$) of $n$:

$\sum _{k=0}^{n-1}\exp \left(- \dfrac{2kr \pi i}{n} \right)=n\delta_{r,nl}$,

so we may multiply by $n$, replace $r$ with $nl$ and sum over $l$:

$\ln(f \left( z \right)) =-\sum _{q=1}^{\infty}\sum _{l=1}^{\infty}\left(\dfrac{z}{{\pi}q}\right)^{2nl}l^{-1} $.

Now recognise the sum over $q$:

$\sum _{q=1}^{\infty}\dfrac{1}{q^{2nl}}=\zeta(2nl)$

and so it follows:

$\ln(f \left( z \right)) =-\sum _{l=1}^{\infty}\zeta(2nl)\left(\dfrac{z}{{\pi}}\right)^{2nl}l^{-1} $.

Now take the exponential:

$f \left( z \right) =\exp\left(-\sum _{l=1}^{\infty}\zeta(2nl)z^{2nl}{\pi}^{-2nl}l^{-1} \right)$,

$f \left( z \right) =\prod _{l=1}^{\infty}\exp\left(-\zeta(2nl)z^{2nl}{\pi}^{-2nl}l^{-1} \right)$,

$f \left( z \right) =\prod _{l=1}^{\infty}\sum _{u=0}^{\infty}\left(-\zeta(2nl)z^{2nl}{\pi}^{-2nl}l^{-1} \right)^u(u!)^{-1}$.

We are interested in the coefficient of $z^{2n}$ and this will not receive any contribution from $1<l$ or $1<u$ and thus:

$f \left( z \right) =1-\zeta(2n)z^{2n}{\pi}^{-2n}+O(z^{4n})$,

and consequently:

$g(z)=2(-1)^n(2n)!2^{-2n}f(z)= 2(-1)^n(2n)!2^{-2n}-2(-1)^n(2n)!2^{-2n}\zeta(2n)z^{2n}{\pi}^{-2n}+O(z^{4n})$.

From which it follows that the coefficient of $z^{2n}$ is:

$-2(-1)^n(2n)!\zeta(2n)(2\pi)^{-2n}=B(2n)$.

As you have supplied an alternative expression for the same coefficient your expression must also equal $B(2n)$.

It seems like an interesting expression, do you have an application for it? Can you use it to calculate Fourier coefficients in products of series for example or something like that?