Asymptotic expansion of $\int_0^{+\infty} \frac{ne^{-\sqrt{t}}}{1+n^2t^2}\,dt$
Changing $u=z^2/n$ and $\zeta=n^{-1/2}$, the integral can be expressed as \begin{align} I(\zeta) &= \int_0^{+\infty} \frac{ne^{-\sqrt{t}}}{1+n^2t^2} \,dt\\ &=\int_0^{+\infty} \frac{2u}{1+u^4}e^{-\zeta u}\,du \end{align} We use the Mellin transform method for Laplace transforms with small parameters (DLMF). Here, \begin{align} h(u)&=\frac{2u}{1+u^4}\\ &\sim 2\sum_{s=0}^\infty (-1)^s u^{-3-4s} \end{align} for $u\to\infty$ and $h(u)=O\left( u \right)$ for $u\to 0$. The Mellin transform \begin{equation} \mathfrak M\left[ h(u)\right](z)=\frac{\pi}{2\sin\left( \frac{\pi}{4}(z+1) \right)} \end{equation} for $-1<z<3$. The asymptotic expansion can then be expressed as \begin{equation} I(\zeta)= \mathfrak M\left[ h(u)\right](1)+\sum_{k=2}^l \operatorname{Res}\left[-\zeta^{z-1}\Gamma(1-z)\frac{\pi}{2\sin\left( \frac{\pi}{4}(z+1) \right)};z=k \right]+O\left( \zeta^{l-\delta-1} \right) \end{equation} for $l>2$ and $\delta$ is a arbitrary small positive constant. Here $\psi(z)=\Gamma'(z)/\Gamma(z)$. We find, the values of the residues $R_k$ \begin{align} R_k&=\frac{2(-1)^{\frac{k+1}{4}}\zeta^{k-1}}{(k-1)!}\left[ \ln \zeta-\psi(k-2)-\frac{2k-3}{(k-1)(k-2)}\right]\text{ for } k= 3+4p \qquad (p=0,1,2\ldots)\\ &=\frac{(-1)^{k-1}\pi\zeta^{k-1}}{2(k-1)!\sin\left( \frac{\pi}{4}(k+1) \right)}\text{ for } k\ge 2 \text{ and } k\ne 3+4p\qquad (p=0,1,2\ldots) \end{align} By keeping the first terms ($k=2,3$) one obtains \begin{equation} I\sim\frac{\pi}{2}-\frac{\pi}{\sqrt{2n} }-\frac{1}{n}\left[ \ln \zeta+\gamma-\frac{3}{2}\right] \end{equation}
What's about good old asymptotic matching? This is at least feasible to get the first two (or maybe three,four,... terms if you are bit more patient) of the expansion, without the need to have heavy artillery from complex analysis at hand.
Let us first fix a constant $\Delta<<1$ depending on $n$ such that $\Delta n \rightarrow \infty $ as $n\rightarrow\infty$. The exact size depends on the order of matching and can be specified in the end so that the Expansion is consistent ($1/n^{\delta},\,\,1>\delta>0$ will suit our needs) Let us write ($J=I/n$)
$$ J=J_1+J_2=\int_0^{\Delta}f_n(x)dx+\int_{\Delta}^{\infty}f_n(x)dx $$
where $f_n(x):=\frac{e^{-\sqrt{x}}}{1+n^2 x^2}$.
Let us start with the Analysis of $J_1$
In the interval $[0,\Delta]$ we can expand the exponential as a Taylor series since $x<<1$:
$$ J_1\sim\int_0^{\Delta}\frac{\color{\red}{1}-\color{\green}{\sqrt{x}}+\color{\orange}{O(x)}}{1+n^2 x^2}dx $$
Integrating termwise yields
$$ J_1\sim\color{\red}{\frac{\arctan(\Delta n)}{n}}+\color{\green}{\frac{e^{i\pi/4}(\arctan(e^{i\pi/4}\sqrt{\Delta n})-\text{arctanh}(e^{i\pi/4}\sqrt{\Delta n}))}{n^{3/2}}}+\color{\orange}{O(\log(\Delta n)/n^2)} $$
Expanding up to $O(\sqrt{\Delta} n^{2})^{-1}$ we get (discarding small terms) $$ J_1\sim\color{\red}{\frac{\pi}{2 n}-\frac{1}{\Delta n^2}}-\color{\green}{\frac{\pi}{\sqrt{2}n^{3/2}}+2\frac{1}{n^2\sqrt{\Delta}}}+\color{\orange}{O(\log(n)n^{-2})} $$
Now, let's consider $J_2$
For $J_2$ we can Taylor the fraction around infinity since $xn>>1$ on $[\Delta,\infty]$
$$ J_2\sim\frac{1}{n^2}\int_{\Delta}^{\infty}\frac{e^{-\sqrt{x}}}{x^2}\left(\color{\purple}{1}+O((nx)^{-2})\right) $$
Integrating by parts twice we get
$$ J_2\sim\color{\purple}{\frac{e^{-\sqrt{\Delta}}}{\Delta n^2}-\frac{e^{-\sqrt{\Delta}}}{\sqrt{\Delta} n^2}} $$
Expanding up to $O(\sqrt{\Delta} n^{2})^{-1}$ yields (discarding small terms) $$ J_2\sim\color{\purple}{\frac{1}{\Delta n^2}-2\frac{1}{n^2\sqrt{\Delta}}} $$
Adding both integrals, we see that the parts containing $\Delta$ exactly cancel which makes our Expansion consistent up to the order in question, and we can state our final result
$$ I=nJ=n(J_1+J_2)\sim \color{\red}{\frac{\pi}{2 }}-\color{\green}{\frac{\pi}{\sqrt{2}n^{1/2}} }+\color{\orange}{O(\log(n)n^{-1})} $$
PS:Thanks for the idiotic downvote
Using a CAS. $$I = \int_0^{+\infty} \frac{ne^{-\sqrt{t}}}{1+n^2t^2} \,dt=\frac{64 \sqrt{2} }{\pi ^{3/2}}n^2\, G_{1,5}^{5,1}\left(\frac{1}{256 n^2}| \begin{array}{c} \frac{3}{2} \\ 1,\frac{5}{4},\frac{3}{2},\frac{3}{2},\frac{7}{4} \end{array} \right)$$ where appears the Meijer G function.
Expanding for large values of $n$ $$I=\frac{\pi }{2}-\frac{\pi }{\sqrt{2}\,n^{1/2}}+\frac{\log \left({n}\right)-2 \gamma +3}{2\, n}+\frac{\pi }{6 \sqrt{2}\,n^{3/2}}-\frac{\pi }{48\, n^2}+O\left(\frac{1}{n^{5/2}}\right)$$ For $n=10$, the "exact" result would be $\approx 1.086826$ while the above gives $\approx 1.086776$.
No problem to get more terms (if you want).