A more general form is known : $\tag{1}f(c)=\sum_{n,m}'\frac 1{(n^2+m^2)^c}=4\,\zeta(c)\,L_{-4}(c)$ with $\ L_{\sigma}(c):=\sum_{n\ge 1} \,\left(\frac {\sigma}n\right)\,n^{-c}\quad$ (where $\ \displaystyle\left(\frac {\sigma}n\right)$ is the Kronecker symbol and the sum is on all signed integers $n,m$ excluding $n=m=0$ as indicated by ' )
For $\;\sigma=-4\ $ only the odd terms remain with alternating signs : $\tag{2}L_{-4}(c)=\sum_{n\ge 1} \,\left(\frac {-4}n\right)\,n^{-c}=\sum_{n\ge 1} \frac {(-1)^{n-1}}{(2n-1)^c}=\beta(c)$ with $\beta$ the Dirichlet beta function so that (distinguishing the values $n=0$ and $m=0$) : $4\,\zeta(c)\,\beta(c)=\sum_{n,m}'\frac 1{(n^2+m^2)^c}=2\sum_{m\ge 1}\frac 1{(0^2+m^2)^c}+2\sum_{n\ge 1}\frac 1{(n^2+0^2)^c}+4\sum_{n,m\ge 1}\frac 1{(n^2+m^2)^c}$ and (dividing by $4$) a generalization of your nice result : $\tag{3}\sum_{n,m\ge 1}\frac 1{(n^2+m^2)^c}=\,\zeta(c)\,\beta(c)-\zeta(2c)$
The (rather non-trivial) derivation of a further generalization of $(1)$ is provided in Borwein, Bailey and Girgensohn's "Experimentation in Mathematics" (page $167-169$ of my $2004$ edition).
They define :
$\tag{4}\zeta(a,b,c) = \sum_{n,m}' \frac{n^{2a}\,m^{2b}}{(n^2+m^2)^c}$
($a,\;b,\;c\;$ are nonnegative integers)
(the sum was for $n,m\ge 0$ but I think that "$\ge 0$" is a typo that doesn't appear page $169$)
The page $168$ of the book contains (I hope that a one page re-transcription will not be a problem, extracts of the book may be viewed at Amazon) :
"The following identity expresses $\zeta(a,0,c)$ by using a Bessel function expansion of the normalized Mellin transform $M_c(f) = \frac 1{\Gamma(c)} \int_0^\infty f(t)\,t^{c-1}\,dt$ of the function $t \to \sum_{n,m} n^{2a}q^{n^2+m^2} =\sum_n n^{2a}q^{n^2} \theta_3(q)$ with $q=e^{-t}$ after using the theta transform $(2.15)$ $\theta_3\bigl(e^{\pi t}\bigr)= \sqrt{\frac{\pi}t}\,\theta_3\left(e^{\frac{\pi}t}\right)$ This leads to the identity $\tag{5}\zeta(a,0,c)= 2\delta_{0a}\zeta(2c)+ 2\beta\left(c-\frac12,\frac12\right) \zeta(2c-2a-1)+4\sum_{p\ge 1} \sigma_{[2a+1-2c]}(p)\,{\cal E}_c(p)$ This is valid for real $ac$ with $d=c-a>1$ and presumably provides an analytic continuation of $\zeta(a,0,c)$ for the sum over positive integers. Here $\sigma$ is a divisor function $\sigma_{[d]}(p) =\sum_{n|p} n^d,$ and ${\cal E}_c(p)= \frac{\sqrt{\pi}}{\Gamma(c)} 2(\pi p)^{c-1/2}K_{\bigl(c-\frac 12\bigr)}(2\pi p)$ is derived from the modified Bessel function of half integer order, $K_{\bigl(c-\frac 12\bigr)}$, of the second kind.
When $c=N$ is integer, ${\cal E}_N(p)$ is of the form $\pi\,e^{-2\pi p} P_N(\pi p)$ where $P_N$ is a rational polynomial of degree $N-1$ with positive coefficients : $P_1(x)=1,\ P_2(x)=x+1/2,\ P_3(x)=x^2/2+3/4x+3/8,\ P_4(x)=1/6x^3+1/2x^2+5/8x+5/16$.
In general $P_N(x)=\sum _{k=0}^{N-1}\frac{{{N+k-1}\choose{N-1}}x^{N-1-k}}{(N-k-1)!4^k}$ This allows one to very efficiently compute $\zeta(a,0,c)$ via $(5)$, using roughly $D/4$ terms for $D$ digits. Note also that for fixed $c$ and variable $a$ only the powers in $\sigma_{[2a+1-2c]}$ vary so most of the computation can be saved.
Then the general integer case follows from $\tag{6}\zeta(a,b,c)=\sum_{k=0}^e (-1)^{e-k}{{e}\choose{k}}\zeta(a+b-k,0,c-k)$ where $\ e=\min(a,b)$."
The formula at the beginning was $\zeta(0,0,c)$ and you asked merely for $\zeta(0,0,2)$. I think there is a much simpler answer for the second one (so that the problem remains...) but generalizations are of value too!
Let's add some keywords and references I just found (MathWorld and others) :
Hoping this helped too,