One possibility is the following:

  1. Replace $$\frac{1}{k^2+a^2}=\int_0^{\infty}e^{-s(k^2+a^2)}ds.$$

  2. The integral over $n$ components of $k$ now factorizes into the product of $n$ Gaussian integrals $$\prod_{\alpha=1}^n\int_{-\infty}^{\infty}e^{-sk_{\alpha}^2+ik_{\alpha}x_{\alpha}}dk_{\alpha}=\prod_{\alpha=1}^n\sqrt{\frac{\pi}{s}}\,e^{-x_{\alpha}^2/4s}=\left(\frac{\pi}{s}\right)^{\frac{n}{2}}e^{-r^2/4s}.$$

  3. Therefore the $n$-fold integral $I(x)$ can be rewritten as an ordinary integral $$I(x)=\int_{0}^{\infty}\left(4\pi s\right)^{-n/2}e^{-a^2 s-r^2/4s}ds,$$ where $r^2=x\cdot x$. Making the change of variables $t=\frac{2a s}{r}$, this can be transformed into $$I(x)=\frac{1}{4\pi}\left(\frac{2\pi r}{a}\right)^{1-\frac{n}{2}}\int_0^{\infty}s^{-\frac{n}{2}}e^{-\frac{ar}{2}\left(s+s^{-1}\right)}ds.$$

  4. Using integral representation of the Macdonald function, the last integral reduces to $$\boxed{I(x)=\frac{1}{2\pi}\left(\frac{2\pi r}{a}\right)^{1-\frac{n}{2}}K_{\frac{n}{2}-1}\left(ar\right)}$$

Remark: For odd dimension $n=2m+1$, Macdonald function $K_{m-\frac12}(ar)$ can be expressed in terms of elementary functions. In particular, $K_{\frac12}(z)=\sqrt{\frac{\pi}{2z}}\, e^{-z}$, which reproduces the formula for $n=3$ mentioned above.