Integrating by parts twice,
$$ \begin{align} \int_{0}^{1} (\arctan x)^{2} \ dx &= x (\arctan x)^{2} \Big|^{1}_{0} - 2 \int_{0}^{1} \frac{x \arctan x}{1+x^{2}} \ dx \\ &= \frac{\pi^{2}}{16} - 2 \int_{0}^{1} \frac{x \arctan x}{1+x^{2}} \ dx \\ &= \frac{\pi^{2}}{16} - \arctan(x) \ln(1+x^{2}) \Big|^{1}_{0} + \int_{0}^{1} \frac{\ln (1+x^{2})}{1+x^{2}} \ dx \\ &= \frac{\pi^{2}}{16} - \frac{\pi}{4} \ln 2 + \int_{0}^{1} \frac{\ln(1+x^{2})}{1+x^{2}} \ dx \end{align}$$
Let $x = \tan t $.
Then
$$\begin{align}\int_{0}^{1} (\arctan x)^{2} \ dx &=\frac{\pi^{2}}{16} - \frac{\pi}{4} \ln 2 - 2 \int_{0}^{\pi /4} \ln (\cos t) \ dt \\ &= \frac{\pi^{2}}{16} - \frac{\pi}{4} \ln 2 -2 \int_{0}^{\pi /4} \left( \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \cos (2nt) - \ln 2 \right) \ dt \\ &= \frac{\pi^{2}}{16} - \frac{\pi}{4} \ln 2 - 2 \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n}\int_{0}^{\pi /4} \cos (2nt) \ dt + \frac{\pi}{2} \ln 2 \\ &= \frac{\pi^{2}}{16} + \frac{\pi}{4} \ln 2 - \sum_{n=1}^{\infty} \frac{\sin \left(\frac{\pi n}{2} \right)}{n^{2}} \\ &= \frac{\pi^{2}}{16} + \frac{\pi}{4} \ln 2 - \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)^{2}} \\ &= \frac{\pi^{2}}{16} + \frac{\pi}{4} \ln 2 - C \end{align}$$
Fourier series of Log sine and Log cos
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*}
Best Answer
Continue with $$I= -2\int_{0}^{1}\frac{\arctan x\ln x}{1+x^{2}}\overset{x\to \frac1x}{dx}= \frac\pi2 \int_1^\infty \frac{\ln x}{1+x^2}dx -\int_0^\infty \frac{\arctan x\ln x}{1+x^2} dx$$ where $\int_1^\infty \frac{\ln x}{1+x^2}dx=G$ and \begin{align} \int_0^\infty \frac{\arctan x\ln x}{1+x^2}dx =& \int_0^\infty \int_0^1 \frac{x\ln x}{(1+x^2)(1+y^2x^2)} \overset{x\to \frac1{xy}}{dx}dy\\ = & \ \frac1{2}\int_0^1\int_0^\infty \frac{-x\ln y}{(1+x^2)(1+{y^2}x^2)} {dx}\ dy\\ =& \ \frac12\int_0^1\frac{\ln^2 y}{1-y^2}dy =\frac78\zeta(3) \end{align} Plug back into $I$ to ontain $$I= \frac\pi2G- \frac78\zeta(3)$$