The following is just a sketch, for the detail you can find in the reference of Zagier [http://people.mpim-bonn.mpg.de/zagier/files/scanned/ValeursZeta/ZetaFunctionRQF.pdf]. Let us denote
$$f(x,y)=(x+y)e^{-(axy+(x-y)^2)}.$$
Then you can check that
$$F(s):=\sum_{m,n\ge 1}\frac{m+n}{(amn+(m-n)^2)^s}=\frac{1}{\Gamma(s)}\int_{0}^{\infty}u^{s-3/2}\left(\sum_{m,n\ge 1}f(\sqrt{u}m,\sqrt{u}n)\right)\,du$$
by direct calculation. Then, from a result of Zagier [page 10th of http://people.mpim-bonn.mpg.de/zagier/files/scanned/ValeursZeta/ZetaFunctionRQF.pdf], actually the Euler–Maclaurin formula, we have the following asymptotic expansion:
\begin{align}
\sum_{m,n\ge 1}f(mt,nt)\sim&\frac{1}{t^2}\int_{{\mathbb{R}}_+^2}f(x,y)\,dx\,dy+\sum_{r,s\ge 0}\beta_r\beta_sf^{(r,s)}(0,0)t^{r+s}\\
&+\frac{1}{t}\sum_{r\ge 0}\beta_rt^r\left(\int_{\mathbb{R}_+}f^{(0,r)}(x,0)\,dx+\int_{\mathbb{R}_+}f^{(r,0)}(0,y)\,dy\right)
\end{align}
for $t\rightarrow 0^+$, where $\beta_r=(-1)^rB_{r+1}/(r+1)!, r\in\mathbb{Z}_{\ge 0}$ and $B_r$ denote the $r$-th Bernoulli number. There we have as $|u|\le 1,$
$$\sum_{m,n\ge 1}f(\sqrt{u}m,\sqrt{u}n)=\frac{I_f}{u}+\frac{c_{-1}}{u^{1/2}}+c_0+O(u^{1/2})$$
with $I_f,c_{-1}, c_0$ be defined as above asymptotic expansion.
Hence we can obtain that
$$F(s)=\frac{1}{\Gamma(s)}\int_{0}^{1}u^{s-3/2}\left(\frac{I_f}{u}+\frac{c_{-1}}{u^{1/2}}+c_0\right)\,du+h(u)$$
with $h(u)$ is an analytic function on $\Re(s)>0$. Moreover, for $\Re(s)>3/2$,
$$F(s)=\frac{1}{\Gamma(s)}\left(\frac{I_f}{s-3/2}+\frac{c_{-1}}{s-1}+\frac{c_0}{s-1/2}\right)+h(s).$$
This means that $F(s)$ can analytic continuation to all $\Re(s)>0$. The following is easy and I omit the detail.