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?
Best Answer
The product for $J(s)$ converges if $\Re s>1$ [@jjagmath is a bit off], and the Jacobi triple product gives \begin{align*} J(s)&=\prod_{p\in\mathbb{P}}\sum_{k\in\mathbb{Z}}p^{-k^2 s} \\&=\prod_{p\in\mathbb{P}}\prod_{n\geqslant 1}(1-p^{-2ns})(1+p^{-(2n-1)s})^2 \\&=\prod_{n\geqslant 1}\prod_{p\in\mathbb{P}}\frac{(1-p^{-2ns})(1-p^{-(4n-2)s})^2}{(1-p^{-(2n-1)s})^2} \\&=\prod_{n\geqslant 1}\frac{\zeta^2\big((2n-1)s\big)}{\zeta(2ns)\zeta^2\big((4n-2)s\big)}. \end{align*} With $\zeta$ analytically continued, the last product converges locally uniformly for $\Re s>0$ and $s$ not a singularity of any term, thus it defines the analytic continuation of $J(s)$ for these values of $s$. The poles of $J(s)$ are nontrivial zeros of $\zeta(ks)$ for even $k$, and we know that any neighborhood of any $s$ with $\Re s=0$ contains infinitely many of the poles. Thus, $\Re s=0$ is the natural boundary of $J(s)$, and the notation "$J(-s)$" makes no sense for $\Re s\geqslant 0$.