For $r$ an integer $\geq 2$ you could start from the closed form expression for $r=2$ and then obtain larger integer values of $r$ by differentiation with respect to $x_0$. The result for $r=2$ isFor $r$ an integer $\geq 2$ you could start from the closed form expression for $r=2$ and then obtain larger integer values of $r$ by differentiation with respect to $x_0$. The result for $r=2$ is
$$ I_2=\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^2} dy=\frac{\sqrt{2 \pi }}{2x_0^2}+\frac{\pi e^{x_0^2/2}}{4x_0^3 e^{y_0^2/2}}(w+\bar{w}), $$ with $z=x_0+iy_0$, $w=(1-x_0z)\text{erfc}(z/\sqrt{2})\exp(i x_0 y_0)$.$$ I_2=\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^2} dy=\frac{\sqrt{2 \pi }}{2x_0^2}+\frac{\pi e^{x_0^2/2}}{4x_0^3 e^{y_0^2/2}}(w+\bar{w}), $$
with $z=x_0+iy_0$, $w=(1-x_0z)\text{erfc}(z/\sqrt{2})\exp(i x_0 y_0)$.
For the large-$x_0$ asymptotics I follow the suggestion of Venkataramana:
$$I_r=\int_{-\infty}^\infty \frac{e^{-y^2/2}}{((y+y_0)^2+x_0^2)^r} dy=\frac{\sqrt{2\pi}}{\Gamma(r)}\int_0^\infty\frac{u^{r-1} }{\sqrt{2 u+1} }e^{-u \left(\frac{y_0^2}{2 u+1}+x_0^2\right)}\,du$$ Now for $x_0\rightarrow\infty$ the integrand contributes mainly near $u=0$, so I may replace $2u+1\mapsto 1$, and then find $$I_r\rightarrow\sqrt{2 \pi } \left(x_0^2+y_0^2\right)^{-r}.$$ Here are some numerical checks, blue is the exact integral, gold the large-$x_0$ form, all plotted as a function of $x_0$ for different values of $y_0$ and $r$.
top row from left to right: $y_0=1,r=2$; $y_0=1,r=4$ bottom row from left to right: $y_0=10,r=2$; $y_0=10,r=4$.