This series is a good candidate for the Poisson summation formula. Let $f(x) = \frac{1}{(1+a^2 x^2)^{3/2}}$ and $f_{\frac{1}{2}}(x) = f(x+1/2)$. Define the sum $J = \sum_{n\in \mathbb{Z}} f_{\frac{1}{2}}(n)$. It is related to the desired sum $I = \sum_{n=0}^{\infty} f_{\frac{1}{2}}(n)$ by the formula $I=J/2$, because $f_{\frac{1}{2}}(-n-1) = f_{\frac{1}{2}}(n)$.
By the Poisson summation formula, we have the equality
$$
J =
\sum_{n\in\mathbb{Z}} f_{\frac{1}{2}}(n) =
\sum_{p\in\mathbb{Z}} F_{\frac{1}{2}}(2\pi p) ,
$$
where $F(k) = \int_{-\infty}^{\infty} dx\, e^{-ikx} f(x)$ and $F_\frac{1}{2}(k) = \int_{\infty}^{\infty} dx\, e^{-ikx} f_{\frac{1}{2}}(x) = e^{ik/2} F(k)$. The idea is that for each larger value of $p$, the terms of the sum in Fourier space will be of larger and larger subleading order in $a$, so that only the first few terms of that serious would be needed to get a good asymptotic estimate.
Now, we can use some integration by parts and complex contour integration methods to evaluate the actual Fourier transform $F(k)$.
\begin{align}
F(k)
&= \int_{-\infty}^{\infty} dx\, e^{-ikx} \frac{1}{(1+a^2 x^2)^{3/2}} \\
&= \int_{-\infty}^{\infty} dx\, e^{-ikx} \frac{d}{dx}\frac{x}{(1+a^2 x^2)^{1/2}} \\
&= ik \int_{-\infty}^{\infty} dx\, e^{-ikx} \frac{x}{(1+a^2 x^2)^{1/2}} \\
\end{align}
The integration by parts helps make the subsequent deformed contour integrals finite. For $k>0$ the integral can be deformed into the lower half complex plane, while for $k<0$ it can be deformed into the upper half. In each case, the new contour will run along both sides of a branch cut, extending from either $x=i/a$ or $x=-i/a$ outward to infinity. Since $f(x)$ is real and symmetric, so is $F(k)$. So it is enough to consider only $k>0$ and the lower branch cut. Note that the imaginary part of the square root will be negative along the branch of that contour running to infinity.
Letting $x=-i(z+1)/a$ and summing the contributions to the contour from either side of the branch cut we get
\begin{align}
F(k)
&= (ik) (-2i/a) e^{-k/a} \int_0^\infty dz\, e^{-kz/a} \frac{-i(z+1)/a}{-i\sqrt{(z+1)^2-1}} \\
&= \frac{2k}{a^2} e^{-k/a} \int_0^\infty dz\, \frac{e^{-kz/a}}{\sqrt{z}} \frac{1+z}{\sqrt{2+z}} \\
&\sim \frac{2k}{a^2} e^{-k/a} \int_0^\infty dz\, \frac{e^{-kz/a}}{\sqrt{z}} \frac{1}{\sqrt{2}} \left(1+ \frac{3}{4}z - \frac{5}{32} z^2 + \cdots \right) \\
&\sim \frac{\sqrt{2\pi}}{a} e^{-k/a} \frac{k^{1/2}}{a^{1/2}} \left(1 + \frac{3}{8}\frac{a}{k} - \frac{15}{128}\frac{a^2}{k^2} + \cdots \right),
\end{align}
where the last two steps constitute an asymptotic series in $a/k$ obtained by expanding the integrand in a non-everywhere uniformly converging power series. Note that the singularity at $z=0$ is of integrable type. It would have been non-integrable without the integration by parts done before deforming the contour. Note also that the above formulas don't work for $k=0$. Fortunately, we can obtain that value directly,
$$
F(0) = \int_{-\infty}^{\infty} dx\, \frac{1}{(1+a^2x^2)^{3/2}}
= \int_{-\infty}^{\infty} dx\, \frac{d}{dx} \frac{x}{(1+a^2x^2)^{1/2}}
= \frac{2}{a} .
$$
So, finally, an $a\sim 0$ asymptotic formula for the desired sum is
\begin{align}
I
&= \frac{1}{2}\left(F(0) + 2\cos(\pi)F(2\pi) + \cdots \right) \\
&\sim \frac{1}{a} - \frac{2\pi}{a^{3/2}} e^{-2\pi/a} + \cdots ,
\end{align}
where more terms could be computed as needed.
