This reminds me of an old AoPS thread. Not exactly the same problem - that was $\sum_k \frac1{x_k^2}$ - but similar ideas should apply. Following my argument over there, we will try to find an analytic function with residues $\frac1{1+x_k^2}=\cos^2 x_k$ at its poles, which also decays in the right places (away from the real axis, and on the vertical segments we'll use).
$$\text{Let }f(z) = \frac{z\sin z}{(z^2+1)(\sin z-z\cos z)}$$
This $f$ has poles at each $x_k$, at $0$, and at $\pm i$. Each of those poles are simple; this is clear for the $x_k$ and $\pm i$, and at zero $f(z)=\frac{z^2+O(z^4)}{(1+O(z^2))(z-\frac16z^3+O(z^5)-z+\frac12z^3+O(z^5))}=\frac{z^2+O(z^4)}{\frac13z^3+O(z^5)}=\frac3z+O(z)$.

Now, the residues. At a root $x_k$, we have $\sin x_k=x_k\cos x_k$, or $x_k=\tan x_k$. The residue of $f$ there is
$$R(x_k) = \left.\frac{z\sin z}{z^2+1}\right|_{z=x_k}\cdot\left.\frac1{\frac{d}{dz}(\sin z-z\cos z)}\right|_{z=x_k} = \frac{x_k\sin x_k}{x_k^2+1}\cdot\frac1{x_k\sin x_k}=\frac1{x_k^2+1}$$
The residue $R(0)$ at zero is $3$ by our power series calculation above. The residues at $\pm i$ are
$$R(\pm i) = \left.\frac{z\sin z}{\sin z-z\cos z}\right|_{z=\pm i}\cdot\left.\frac1{\frac{d}{dz}(z^2+1)}\right|_{z=\pm i} = \frac{\pm i}{1-\pm i\cot\pm i}\cdot \frac1{\pm 2i}=\frac1{2(1-\coth 1)}$$
using $\cos it=\cosh t$ and $\sin it=i\sinh t$ there.

And now, we build a contour integral to take advantage. Consider $I(N,M)=\int_C f(z)\,dz$, where $f$ is the rectangle with vertices at $N\pi+iM, -N\pi+iM, -N\pi-iM, N\pi-iM$ for some large $M$ and large integer $N$. For convenience, also rewrite $f(z)=\frac{z}{z^2+1}\cdot \frac1{1-z\cot z}$. Now, on a vertical segment $N\pi+iy$, $\cot z=\frac{\cos(N\pi)\cosh y}{i\cos(N\pi)\sinh y}=-i\coth y$ and
$$f(N\pi+iy)= \frac{N\pi + iy}{N^2\pi^2-y^2+1+2N\pi yi}\cdot \frac1{N\pi i+1-y\coth y}$$
As long as $M \le N\pi$, $N^2\pi^2>y^2$ and $|z^2+1|>z^2$ on the vertical segment. The first term is then less than $|z|^{-1}<\frac1{N\pi}$ in absolute value, and the second term is less than $\frac1{N\pi}$ in absolute value. Integrating that product over a length of $2M$, we get at most $\frac{2M}{(N\pi)^2}\le\frac{2}{N\pi}$.

On the horizontal segments, we have that $\cot (x+iy)$ tends to $-i$ uniformly as $y\to\infty$ and to $i$ uniformly as $y\to-\infty$. Choose $M$ so that it's within $\epsilon$ on the segments $x+\pm iM$, and we have
$$|f(z)|\le \frac{|z|}{|z|^2-1}\cdot\frac1{(1-\epsilon)|z|-1}$$
By limit comparison to $\frac1{M^2}$, the integral of this is at most a constant times $\frac{N}{M^2}$. As long as $N$ isn't too much larger than $M$ - say, $M=N\pi$ to match the other direction - this goes to zero.

Therefore, $\lim_{N\to\infty}I(N,N\pi) = 0$. And now, we calculate the integral by the residue theorem. Inside the contour, there are $2N+1$ poles: $0,i,-i,x_1,-x_1,\dots,x_{N-1},-x_{N-1}$. The residue at each $x_k$ is equal to that at each $-x_k$, and
$$I(N,N\pi) = 2\pi i\left(3+\frac1{2(1-\coth 1)}+\frac1{2(1-\coth 1)} + \sum_{k=1}^{N-1}\frac1{x_k^2+1}+\sum_{k=1}^{N-1}\frac1{x_k^2+1}\right)$$
$$0 = 3+\frac1{1-\coth 1}+2\sum_{k=1}^{\infty}\frac1{x_k^2+1}$$
$$\sum_{k=1}^{\infty}\frac1{x_k^2+1} = -\frac12\left(\frac{\sinh 1}{\sinh 1-\cosh 1}+3\right)=\frac12\left(e\cdot\frac{e-e^{-1}}{2}-3\right)=\frac{e^2}{4}-\frac74$$