This math.SE answer calculated the following (inverse) Fourier transform in $n$-dimensions
$$ f(x) = \int_{\mathbb{R}^n} d^nk \frac{e^{ik\cdot x}}{|k|^2 + 1} \propto |x|^{-\frac{n}{2}+1} K_{\frac{n}{2}-1}(|x|) \tag{1} $$
where $K_\alpha(x)$ is the modified Bessel function of the second kind. But in some physics applications, one is only interested in the large-$|x|$ asymptotic behavior of $f(x)$. Using the asymptotic expansion
$$ K_\alpha(x) \propto x^{-1/2} e^{-x} [1 + O(x^{-1})] \quad \text{when } \mathrm{arg}(x) < 3\pi/2 $$
we obtain (which reproduces Eq. (14.23) in Assa Auerbach's Interacting Electrons and Quantum Magnetism)
$$ f(x) \sim |x|^{-(n-1)/2} e^{-x} \tag{2} $$
My question is: is there a simpler way than the cited math.SE answer that directly aims to find the asymptotic behavior of $f(x)$ for large $|x|$?
$$f(x)=\int_{-\infty}^\infty...\int_{-\infty}^\infty\frac {e^{ik_1x_1+ik_2x_2+...+ik_nx_n}}{1+k_1^2+...+k_n^2}dk_1...dk_n$$ $$=\int_{-\infty}^\infty...\int_{-\infty}^\infty e^{ik_1x_1+ik_2x_2+...+ik_nx_n}dk_1...dk_n\int_0^\infty e^{-t(1+k_1^2+...+k_n^2)}dt$$ $$\int_0^\infty e^{-t}dt\int_{-\infty}^\infty e^{-t\big((k_1^2-\frac{ix_1}{2t})^2+\frac{x_1^2}{4t^2}\big)}dk_1...\int_{-\infty}^\infty e^{-t\big((k_n^2-\frac{ix_n}{2t})^2+\frac{x_n^2}{4t^2}\big)}dk_n$$ Denoting $x=\sqrt{x_1^2+...+x_n^2}$ and integrating with respect $k_1\,...\,k_n$ $$=\pi^\frac n2\int_0^\infty e^{-t-\frac{x^2}{4t}}t^{-\frac n2}dt\overset{t=\frac{xs}2}{=}\pi^\frac n2\left(\frac x2\right)^{1-\frac n2}\int_0^\infty e^{-\frac x2(s+\frac1s)}s^{-\frac n2}ds$$ Using Laplace' method we can find as many asymptotic terms as we want. The extremum point is $s=1$. For the leading term we get $$f(x)\sim\pi^\frac n2\left(\frac x2\right)^{1-\frac n2}e^{-x}\int_{-\infty}^\infty e^{-\frac x2(s-1)^2}ds=\pi^\frac {n+1}2\left(\frac x2\right)^{\frac{1-n}2}e^{-x}$$