I think I am starting to understand what is going on on. For any Dirichlet character $\chi$, the Eisenstein series
$$\hat{E}_{k,\chi}(z)=-\frac{B_{k,\chi}}{2k}+\sum_{n=1}^{\infty}\left(\sum_{d|n}\chi(d)d^{k-1}\right)q^n$$
is an eigenform of the space $M_k(N,\chi)$. Thus, letting $\chi=\chi_2$ be the unique Dirichlet character modulo $2$ (i.e $\chi_2(n)=n\,\,\mathrm{mod}\,\,2$) we get that
\begin{align*} \hat{E}_{k+1,\chi_2}(z)&=-\frac{B_{k+1,\chi_2}}{2(k+1)}+\sum_{n=1}^{\infty}\left(\sum_{d|n}\chi_2(d)d^{k}\right)q^n\\ &=-\frac{B_{k+1,\chi_2}}{2(k+1)}+\sum_{n=1}^{\infty}\sigma^{(o)}_{k}(n)q^n \end{align*}
Thus, we can differentiate our functions $f_k(z)$ repeatedly $k$ times to get that
\begin{align*} f^{(k)}_k(z)&=\frac{d^k}{dz^{k}}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}\sigma_k^{(o)}(n)}{n^k}e^{z\pi i n}\\ &=(i\pi)^k\sum_{n=1}^{\infty}(-1)^{n-1}\sigma_k^{(o)}(n)e^{z\pi i n}\\ &=-(i\pi)^k\sum_{n=1}^{\infty}\sigma_k^{(o)}(n)e^{(z-1)\pi i n}\\ \end{align*}
and thus
$$\frac{-f^{(k)}_k(z+1)}{(i\pi)^k}=\hat{E}_{k+1,\chi_2}(z)+\frac{B_{k+1,\chi_2}}{2(k+1)}$$
In a way this explains the connection we are seeing, but in another it makes it even more mysterious since there is no reason that the iterated integral of an eigenform should retain its nice properties.
Another (possibly more enlightening) way to think about what is going on is that $f_k(z)$ is the function generated by the $q$-series $\sum_{d|n}\chi_2\left(\frac{n}{d}\right)d^{-k}$, and so perhaps for any choice of Dirichlet character $\chi$ the $q$-series generated by $\sum_{d|n}\chi_2\left(\frac{n}{d}\right)d^{-k}$ will have modular form-like properties, in the sense that
$$z^{\frac{1-k}{2}}g(z)+z^{\frac{k-1}{2}}g\left(\frac{1}{z}\right)=\mathrm{poly}(z)$$
where $\mathrm{poly}(z)$ is a polynomial of degree at most $\frac{k+1}{2}$.