I tried 2 ways to find a closed form, though unsuccessful up to this point.
1st trial. Let $I$ denote the integral, and write
$$ I = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{e^{-x}(1 - e^{-(2n+1)x})}{x (1 + e^{-x})^{2}} \, dx. $$
In order to evaluate the integral inside the summation, we introduce new functions $I(s)$ and $J_n(s)$ by
$$ I(s) = 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \int_{0}^{\infty} \frac{x^{s-1} e^{-x}(1 - e^{-(2n+1)x})}{(1 + e^{-x})^{2}} \, dx =: 8 \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} J_n(s) $$
so that $I = I(0)$. Then it is easy to calculate that for $\Re(s) > 1$, $J_n(s)$ is written as
$$ J_n(s) = \Gamma(s) \left( \eta(s-1) + \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} - (2n+1) \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s}} \right), $$
where $\eta$ is the Dirichlet eta function. Plugging this back to $I(s)$ and manipulating a little bit, we obtain
$$ I(s) = 8\Gamma(s) \left( \frac{\pi}{4} \eta(s-1) - 4^{-s}\left( \zeta(s, \tfrac{1}{4}) - \zeta(s, \tfrac{1}{2}) \right) + \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1} \sum_{k=2n+1}^{\infty} \frac{(-1)^{k-1}}{k^{s-1}} \right). $$
This is valid for $\Re(s) > 1$. But if we can somehow manage to find an analytic continuation of the last summation part, then we may find the value of $I = I(0)$.
2nd trial. I began with the following representation
\begin{align*}
I
&= -2 \int_{0}^{\infty} \frac{1-e^{-t}}{1+e^{-t}} \left( \frac{1}{\cosh t} - \frac{2}{t} ( \arctan(1) - \arctan (e^{-t})) \right) \, \frac{dt}{t} \\
&= \sum_{n=0}^{\infty} \frac{(-1)^{n}}{(2n+1)(2n+2)} \int_{-\infty}^{\infty} \frac{\tanh^{2(n+1)} x}{x^{2}} \, dx.
\end{align*}
With some residue calculation, we can find that
\begin{align*}
\int_{-\infty}^{\infty} \frac{\tanh^{2n} x}{x^{2}} \, dx
&= \frac{2}{i\pi} \, \underset{z=0}{\mathrm{Res}} \left[ \psi_{1}\left(\tfrac{1}{2} + \tfrac{1}{i\pi} z\right) \coth^{2n} z \right] \\
&= 2^{2n+3} \sum_{m=1}^{n} (-1)^{m-1}m (1-2^{-2m-1}) A_{n-m}^{(2n)} \, \frac{\zeta(2m+1)}{\pi^{2m}},
\end{align*}
where $A_m^{(n)}$ is defined by the following combinatoric sum
$$ A_m^{(n)} = \sum_{\substack{ k_1 + \cdots + k_n = m \\ k_1, \cdots, k_n \geq 0 }} \frac{B_{2k_1} \cdots B_{2k_n}}{(2k_1)! \cdots (2k_n)!} = 2^{-2m} [z^{2m}](z \coth z)^{n} \in \Bbb{Q}, $$
where $B_k$ are Bernoulli numbers. Still the final output is egregiously complicated, so I stopped here.
3rd trial. The following yet another representation may be helpful, I guess.
$$ I = \int_{0}^{1/2} \frac{1 - \cot(\pi u/2)}{2} \left\{ \psi_1\left(\tfrac{1+u}{2}\right) - \psi_1\left(\tfrac{1-u}{2}\right) \right\} \, du. $$
We can use these identities:
\begin{align*}
&\mathrm{Li}_2(x)+\mathrm{Li}_2(-x)=\frac{1}{2} \mathrm{Li}_2(x^2) \tag{1} \\
&\mathrm{Li}_2(x)+\mathrm{Li}_2(1-x)=\frac{\pi^2}{6}-\log(x)\log(1-x)\tag{2} \\
&\mathrm{Li}_2(x)+\mathrm{Li}_2\left(\frac{1}{x}\right)=-\frac{\pi^2}{6}-\frac{\log^2(-x)}{2} \tag{3}.
\end{align*}
We can suppose that $x>1$. Additionally, take the principal value of $\log(x)=\ln|x|+i\arg(x)$ i.e. $-\pi<\arg(x)\leq \pi$.
Then equations (3) and (2) become:
$$ \mathrm{Li}_2(x) + \mathrm{Li}_2\left(\frac{1}{x}\right) = \frac{\pi^2}{3}-\frac{1}{2}\ln^2(x)-i\pi\ln(x) \tag{4}$$
$$\mathrm{Li}_2(x)+\mathrm{Li}_2(1-x)=\frac{\pi^2}{6}-\ln(x)\ln(x-1)-i\pi\ln(x) \tag{5}$$
Adding equations (1) and (5):
$$\mathrm{Li}_2(-x)+\mathrm{Li}_2(1-x)=-2\mathrm{Li}_2(x)+\frac{1}{2} \mathrm{Li}_2(x^2)+\frac{\pi^2}{6}-\ln(x)\ln(x-1)-i\pi\ln(x) \tag{6}$$
From equation (4):
$$\mathrm{Li}_2(x) = -\mathrm{Li}_2\left(\frac{1}{x}\right)+\frac{\pi^2}{3}-\frac{1}{2}\ln^2(x)-i\pi\ln(x) \tag{7} $$
$$\mathrm{Li}_2(x^2) = -\mathrm{Li}_2\left(\frac{1}{x^2}\right)+\frac{\pi^2}{3}-2\ln^2(x)-2i\pi\ln(x) \tag{8}$$
Substituting (7) and (8) in (6):
$$ \mathrm{Li}_2(-x) + \mathrm{Li}_2(1-x) = -\frac{1}{2}\mathrm{Li}_2\left(\frac{1}{x^2}\right)+2\mathrm{Li}_{2}\left(\frac{1}{x}\right)-\ln(x)\ln(x-1)-\frac{\pi^2}{3}$$
Hence, adding $\ln(x)\ln(x+1)$:
$$ \mathrm{Li}_2(-x) + \mathrm{Li}_2(1-x) + \ln(x)\ln(x+1) = -\frac{1}{2}\mathrm{Li}_2\left(\frac{1}{x^2}\right)+2\mathrm{Li}_{2}\left(\frac{1}{x}\right)+\ln(x)\ln\left(\frac{1+x}{x-1}\right)-\frac{\pi^2}{3}$$
Given that $\mathrm{Li}_{2}$ is continuous $\forall z\in\mathbb{C}\setminus(1,\infty)$ then $\displaystyle \mathrm{Li}_{2}\left(\frac{1}{x^2}\right),\mathrm{Li}_{2}\left(\frac{1}{x^2}\right) \to 0$ as $x \to \infty$
Additionally, due to the L'Hopital rule:
$$ \lim_{x\to \infty} \ln(x)\ln\left(\frac{1+x}{x-1}\right) = \lim_{x\to \infty} \frac{\ln\left(\frac{1+x}{x-1}\right)}{\frac{1}{\ln(x)}} = 0 $$
$$ \boxed{\lim_{x \to \infty}\left[\mathrm{Li}_2(-x) + \mathrm{Li}_2(1-x) + \ln(x)\ln(x+1) \right]= -\frac{\pi^2}{3}}$$
Best Answer
Not a complete answer, but I'm sure it'll help:
Let $$I(a) = \int_0^\infty e^{-nx}x^a\,dx$$ $$\implies \frac{d^nI(a)}{da^n} = \int_0^\infty e^{-nx}x^a(\ln x)^n\,dx$$ Put $nx \rightarrow v$ in the first integral to get: $$I(a) = \frac1{n^{1+a}}\int_0^\infty e^{-v}v^a\,dv$$ $$\implies I(a) = \frac{\Gamma(1+a)}{n^{1+a}}$$ Now $$\implies \frac{d^nI(a)}{da^n}\bigg|_{a=0} = \frac{d^n}{da^n}\left(\frac{\Gamma(1+a)}{n^{1+a}}\right)\bigg|_{a=0}$$ Which evaluates to: $$\frac1{n}\sum_{k=0}^n(-1)^k\binom{n}{k}\Gamma^{(n-k)}(1+a)\ln^k(n)\bigg|_{a=0}$$ Where $\Gamma^{(n-k)}(1+a)$ is the $(n-k)$th derivative of the Gamma function.