Here is another solution: Let $(I_n)$ by
$$ I_n = \int_{0}^{\infty} \frac{\arctan x}{x(1+x^2)^n} \, dx = \int_{0}^{\frac{\pi}{2}} \frac{\theta}{\sin\theta} \cos^{2n-1}\theta \, d\theta. $$
Then by a simple calculation,
$$ I_n - I_{n+1} = \int_{0}^{\frac{\pi}{2}} \theta \sin\theta \cos^{2n-1}\theta \, d\theta = \frac{1}{2n} \int_{0}^{\frac{\pi}{2}} \cos^{2n}\theta \, d\theta. $$
Since $I_n \to 0$ as $n \to \infty$, we find that
$$ I_n = \sum_{k=n}^{\infty} \frac{1}{2k} \int_{0}^{\frac{\pi}{2}} \cos^{2k}\theta \, d\theta. $$
Splitting the summation as $\sum_{k=n}^{\infty} = \sum_{k=1}^{\infty} - \sum_{k=1}^{n-1}$, we find that
$$ I_n = \frac{\pi}{2}\left( \log 2 - \sum_{k=1}^{n-1} \frac{1}{2k} \frac{(2k-1)!!}{(2k)!!} \right), $$
where $n!!$ denotes the double factorial.
Edit 1. In general, we have
$$ \int_{0}^{\infty} \frac{\arctan^s x}{x(1+x^2)^{n+1}} \, dx = \int_{0}^{\frac{\pi}{2}} \theta^s \cot \theta \, d\theta - \sum_{k=1}^{n} \int_{0}^{\frac{\pi}{2}} \theta^s \sin\theta \cos^{2k-1}\theta \, d\theta. \tag{1} $$
Currently I have no idea how to obtain a simple formula for the following integral
$$ \int_{0}^{\frac{\pi}{2}} \theta^s \sin\theta \cos^{2k-1}\theta \, d\theta, \tag{2} $$
even when $s = 2$. On the other hand, for any $s > 0$ and $N \geq \lfloor s/2 \rfloor$ we have
\begin{align*}
\int_{0}^{\frac{\pi}{2}} \theta^s \cot \theta \, d\theta
&= 2^{-s}\cos\left(\frac{\pi s}{2}\right)\Gamma(1+s)\zeta(1+s) \\
&\quad + \left(\frac{\pi}{2}\right)^s \sum_{k=0}^{N} (-1)^k \pi^{-2k} \frac{\Gamma(2k-s)}{\Gamma(-s)} \eta(2k+1) \\
&\quad + \frac{(-1)^{N+1}}{2^s \Gamma(-s)} \int_{0}^{\infty} \frac{t^{2N+1-s}}{1+t^2} \left( \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n^{1+s}} e^{-\pi n t} \right) \, dt,
\end{align*}
where $\eta(s)$ denotes the Dirichlet eta function. (My solution is somewhat involved, so I will post later if it seems useful to our problem.) In particular, when $s$ is a positive integer, then the integral part vanishes and the formula becomes much simpler. Thus the formula (*) gives a closed form as long as we can figure out the integral (2).
Example 1. For example, when $s = 2$ then we can use $N = 1$ and then
\begin{align*}
\int_{0}^{\frac{\pi}{2}} \theta^2 \cot \theta \, d\theta
&= -\frac{1}{2}\zeta(3) + \frac{\pi^2}{4} \log 2 - \frac{1}{2}\eta(3) \\
&= \frac{\pi^2}{4}\log 2 - \frac{7}{8}\zeta(3).
\end{align*}
Since we can figure out the integral (2) for $s = 2$ and $k = 1, \cdots, 4$, we easily obtain OP's last identity.
Here is another example:
Example 2. Using the formula with $s = 6$, we can check that
\begin{align*}
\int_{0}^{\infty} \frac{\arctan^6 x}{x(1+x^2)^3} \, dx
&= \frac{\pi^6}{64} \log 2 -\frac{45 \pi^4}{128} \zeta(3) + \frac{675 \pi^2}{128} \zeta(5) -\frac{5715}{256} \zeta(7) \\
&\quad - \frac{11 \pi^6}{2048} + \frac{705 \pi^4}{4096} - \frac{8595 \pi^2}{4096} + \frac{135}{16} \\
&\approx 0.0349464822054751922142122595622\cdots.
\end{align*}
For any $a>0$ we have that $f_a(x)=\frac{a}{x^2+a^2}$ is an even function with an absolute maximum at the origin such that $\int_{-\infty}^{+\infty}f_a(x)\,dx=\pi$. In terms of the Laplace transform
$$(\mathcal{L}^{-1} f_a)(s)=\sin(as)$$
holds, and if we assume that
$$ g(x) \stackrel{L^2}{=} c_0 + \sum_{n\geq 1}\left( c_n \cos(nx) + s_n \sin(nx)\right) $$
has a uniformly convergent Fourier series (which is granted, for instance, by $\max(s_n,c_n)=O\left(\frac{1}{n^2}\right)$) we have
$$\begin{eqnarray*} \int_{-\infty}^{+\infty}g(x)f_a(x)\,dx &=& \int_{0}^{+\infty}(g(x)+g(-x))f_a(x)\,dx\\&=&\pi c_0+2\sum_{n\geq 1}c_n\int_{0}^{+\infty}\cos(nx)\frac{a}{a^2+x^2}\,dx\\
&=&\pi c_0+2\sum_{n\geq 1}c_n \int_{0}^{+\infty}\frac{s}{n^2+s^2}\sin(as)\,ds
\\&=&\pi c_0+2\sum_{n\geq 1}c_n \int_{0}^{+\infty}\frac{s}{1+s^2}\sin(nas)\,ds\\&=&\pi c_0+\pi\sum_{n\geq 1}c_n e^{-na}\end{eqnarray*} $$
by the self-adjointness of the Laplace transform. Since $c_n=O\left(\frac{1}{n^2}\right)$, by the dominated convergence theorem
$$ \lim_{a\to 0^+}\int_{-\infty}^{+\infty}g(x)f_a(x)\,dx = \pi\sum_{n\geq 0}c_n = \pi g(0)$$
and by translation
$$ \lim_{a\to 0^+}\int_{-\infty}^{+\infty}g(x-\pi/3)f_a(x-\pi/3)\,dx = \pi g(\pi/3).\tag{1}$$
In order to apply $(1)$, it only remains to prove that $g(x)=\frac{\cos^4(x+\pi/3)}{2+\cos(x+\pi/3)}$ meets the wanted constraints. In order to do that, it is sufficient to estimate the Fourier coefficients of $\frac{1}{2+\cos(x)}$ or $\frac{1}{2-\cos(x)}$. For any $R>1$ we have
$$ \frac{1}{1-\frac{1}{R}e^{ix}}\cdot \frac{1}{1-\frac{1}{R}e^{-ix}}=\frac{R^2}{(1+R^2)-2R\cos(x)}=\sum_{n\geq 0}\frac{e^{nix}}{R^n}\sum_{m\geq 0}\frac{e^{-mix}}{R^m} \tag{2}$$
hence by picking $R=2+\sqrt{3}$ we have
$$\frac{2+\sqrt{3}}{2}\cdot\frac{1}{2-\cos(x)}=\sum_{n\geq 0}\frac{e^{nix}}{(2+\sqrt{3})^n}\sum_{m\geq 0}\frac{e^{-mix}}{(2+\sqrt{3})^m} \tag{3}$$
and the coefficients of the Fourier series of $\frac{1}{2\pm\cos(x)}$ are given by explicit convolutions.
They decay like $\frac{1}{(2+\sqrt{3})^n}$, i.e. way faster than $\frac{1}{n^2}$, and this proves that
$$ \lim_{a\to 0^+}\int_{-\infty}^{+\infty}\frac{\cos^4(x)}{2+\cos(x)}\cdot\frac{a}{(x-\pi/3)^2+a^2}\,dx = \pi\cdot\frac{\cos^4(\pi/3)}{2+\cos(\pi/3)}=\frac{\pi}{40}.\tag{4}$$
Best Answer
Using OP's series expansion, we have
\begin{align*} I &= 2 \sum_{m, n \geq 0} \frac{(-1)^{m+n}}{(2m+1)(2n+1)} \int_{0}^{\infty} \frac{\sin[(2m+1)ax]\sin[(2n+1)ax]}{1+x^2} \, dx \\ &= \pi \sum_{m, n \geq 0} \frac{(-1)^{m+n}}{(2m+1)(2n+1)} \left( e^{-a|2m-2n|} - e^{-a(2m+2n+2)} \right). \end{align*}
We will split the sum into several parts and analyze them separately.
It is straightforward that
$$ \sum_{m, n \geq 0} \frac{(-1)^{m+n}}{(2m+1)(2n+1)} e^{-a(2m+2n+2)} = \arctan^2(e^{-a}). $$
Using the identity $\frac{1}{(2m+1)(2n+1)} = \frac{1}{(2m-2n)(2n+1)} + \frac{1}{(2n-2m)(2m+1)}$, we have
\begin{align*} \sum_{m, n \geq 0} \frac{(-1)^{m+n}}{(2m+1)(2n+1)} e^{-a|2m-2n|} &= \sum_{m = 0} \frac{1}{(2m+1)^2} + \sum_{\substack{m, n \geq 0 \\ m \neq n }} \frac{(-1)^{m-n}}{(m-n)(2n+1)} e^{-a|2m-2n|} \\ &= \frac{\pi^2}{8} + \sum_{n \geq 0} \frac{1}{2n+1} \sum_{\substack{k \geq -n \\ k \neq 0}} \frac{(-1)^k}{k} e^{-2a|k|}. \end{align*}
The last sum can be simplified further: using the fact that $\frac{(-1)^k}{k}e^{-2a|k|}$ is an odd function of $k$,
\begin{align*} \sum_{n \geq 0} \frac{1}{2n+1} \sum_{\substack{k \geq -n \\ k \neq 0}} \frac{(-1)^k}{k} e^{-2a|k|} &= \sum_{n \geq 0} \frac{1}{2n+1} \sum_{k = n+1}^{\infty} \frac{(-1)^k}{k} e^{-2a|k|} \\ &= \sum_{k = 1}^{\infty} \left( \sum_{n=0}^{k-1} \frac{1}{2n+1} \right) \frac{(-1)^k}{k} e^{-2a|k|}. \end{align*}
Symmetrizing the inner sum, we find that
\begin{align*} \sum_{n=0}^{k-1} \frac{1}{2n+1} = \frac{1}{2} \sum_{\substack{i, j \geq 0 \\ i+j = k-1}} \left( \frac{1}{2i+1} + \frac{1}{2j+1} \right) = \sum_{\substack{i, j \geq 0 \\ i+j = k-1}} \frac{k}{(2i+1)(2j+1)}. \end{align*}
Plugging this back,
\begin{align*} \sum_{n \geq 0} \frac{1}{2n+1} \sum_{\substack{k \geq -n \\ k \neq 0}} \frac{(-1)^k}{k} e^{-2a|k|} &= \sum_{i, j \geq 0} \frac{1}{(2i+1)(2j+1)} (-1)^{i+j+1} e^{-a(2i+2j+2)} \\ &= - \arctan^2(e^{-a}). \end{align*}
Combining altogether, we obtain the desired answer.