Roots of some modified Bernoulli polynomials Update
The polynomials are generated as follows: 
Where
$B_n(x) = \sum_{k=0}^n {n \choose k} b_{n-k} x^k$
is used to generate standard Bernoulli polynomials, top plot is generated as follows:
$B_n^{(1)}(x) = -\sum_{k=0}^n {n \choose k} |\Re(b^{(1)}_{n-k})+
\Im(b^{(1)}_{n-k})| x^k$
using definition of $b^{(1)}_{n}$ as outlined below.
(Note: Where $b^{(1)}_{1}={\tilde{\infty}}$, substitute for usual $-\frac{1}{2}$.)
The second plot is exactly the same as above except that the first term of $B_n^{(1)}(x)$ has a positive  value.
Original question
Since
$\zeta(2n) = (-1)^{n+1} \, \frac{2^{2n-1} \, b_{2n} \, \pi^{2n}}{(2n)!},\quad n \ge 0, \,$
where $b_{2n}$ are the Bernoulli numbers, rearranging & substituting $2n$ for $n$, we get:
$b_{n}^{(1)} = \frac{\zeta(n)\, n!}{(-1)^{\frac{n}{2}+1} \, 2^{n-1} \, \pi^{n}}$
Using $b_{n}^{(1)}$ to generate a series similar to the usual Bernoulli polynomials (the main difference being that the odd real $b_{n}$ are give complex values rather than 0), the roots of $B_{n}^{(1)}(x)$ differ from the usual pattern of the $B_{n}(x)$ roots and fit a scaled Szegő Curve even more tightly (bar a couple of spare roots off the diagram) than the roots of the partial sums of $e^x$ (white, outlined points):

Another curious modfication, with roots plotted against scaled Szegő Curve, circle (radius $\frac{n}{e}$) & roots of $e^x - 1$  (white, outlined points):

Image of the two root plots against a dynamic Szegő Curve, as given by Antonio Vargas' code below (included for comparison purposes).

Despite fairly extensive searches on the internet, I have found no literature on the subject. I was wondering if anyone knew of any?
 A: There are some possibilities which could explain the behavior you notice, but I'm just making guesses based on related phenomena.  Proving any of this is likely difficult.
Just as the usual Bernoulli polynomials approximate the sine and cosine (see, e.g., Theorem 1 in this paper), your polynomials may approximate the exponential function.  In fact they may approximate the partial sums of the exponential function in some sense, which could explain the behavior of their zeros.
I gave some references for the zeros of the usual Bernoulli polynomials which may be enlightening in this answer.
If I understand correctly, you are considering the zeros of $B_n^{(1)}(x) + \text{const.}$ in the second plot.  This mirrors the behavior of the zeros of the partial sums of the exponential function $s_n[\exp](z)$.  The zeros of $s_n[\exp](nz)$ converge to the usual Szegő curve, and the zeros of $s_n[\exp](z) + \text{const.}$ converge to the Szegő curve in the right half plane and to the circle $|z| = 1/e$ in the left half plane.  The zeros on the imaginary axis approximate the zeros of $e^{nz} + \text{const.}$
See, for example, this question on MathOverflow.
This is very much related to the fact that $e^{nz} + c \sim e^{nz}$ for $\Re(z) > 0$ and $e^{nz} + c \sim c$ for $\Re(z) < 0$, though I'm not sure what the link would be to your case.

Below is some Mathematica code to plot the intermediate Szegő curves for the exponential function as mentioned in the comments.
s[n_, z_] := Sum[z^k/k!, {k, 0, n}];  
m = 20;  
Show[  
    ContourPlot[  
        Abs[(x+I y) E^(1-x-I y)]^m == Sqrt[2 Pi m] Abs[(1-x-I y)/(x+I y)],  
        {x, -1, 2}, {y, -3/2, 3/2}],  
    ListPlot[
        {Re[z], Im[z]} /. NSolve[s[m, m z] == 0, z],  
        PlotStyle -> Black]  
]

This produces the output

