Here's a fun proof for the second part of your question. (Your conjecture is off by a factor of $1/2$.)
$$
\bbox[15px,#E6FFF3 ,border:5px groove #FF8B00 ]{\int_{0}^{1}\frac{x^{3}\operatorname{artanh}x}{1+x^{4}}dx=\frac{\pi^{2}}{32}-\frac{1}{8}\ln^{2}\left(1+\sqrt{2}\right)}
$$
Proof. Let the integral equal $\mathcal{I}$. Then map $x \mapsto \tanh x$ so that
$$
\begin{align}
\mathcal{I} &:= \int_{0}^{1}\frac{x^{3}\operatorname{artanh}x}{1+x^{4}}dx \\
&= \int_{0}^{\infty}\frac{x\tanh\left(x\right)^{3}}{1+\tanh\left(x\right)^{4}}\operatorname{sech}^{2}xdx \\
&= \int_{0}^{\infty}\frac{2xe^{2x}\left(e^{x}-1\right)^{3}\left(e^{x}+1\right)^{3}}{\left(e^{2x}+1\right)\left(6e^{4x}+e^{8x}+1\right)}dx \\
&= \int_{\mathbb{R}}\frac{xe^{2x}\left(e^{x}-1\right)^{3}\left(e^{x}+1\right)^{3}}{\left(e^{2x}+1\right)\left(6e^{4x}+e^{8x}+1\right)}dx \\
\end{align}
$$
where in the last line, we use the fact that the integrand is an even function.
Let $z \mapsto \displaystyle \frac{z^{2}e^{2z}\left(e^{z}-1\right)^{3}\left(e^{z}+1\right)^{3}}{\left(e^{2z}+1\right)\left(6e^{4z}+e^{8z}+1\right)}$ and define it by the holomorphic function $\displaystyle f: \mathbb{C} \backslash \left\{n \in \mathbb{Z}: z_1(n), z_2(n), z_3(n)\right\} \to \mathbb{C}$. The simple poles excluded from the domain are
$$
\begin{align}
z_{1}(n)&=\frac{i\pi}{2}\left(2n+1\right) \\
z_{2}(n)&=\frac{1}{4}\ln\left(3-2\sqrt{2}\right)+\frac{i\pi}{4}\left(2n+1\right) \\
z_{3}(n)&=\frac{1}{4}\ln\left(3+2\sqrt{2}\right)+\frac{i\pi}{4}\left(2n+1\right)
\end{align}
$$
which are obtained from finding the zeroes of $\displaystyle \frac{1}{f}$. We also construct a rectangular box contour $C$ with vertices $(-R,0)$, $(R,0)$, $(R,R+i\pi)$, and $(-R,-R+i\pi)$ where $R$ is sufficiently large such that $R \gg \displaystyle \frac{1}{4}\ln(3+2\sqrt{2})$. Below, I provide a visual I made of what $C$ looks like in the phase plot of $f$ with domain coloring and shading.
We write each contribution of the integral over $C$ as
$$
\require{cancel}\oint_C f(z)dz=\cancelto{0}{\int_{-R}^{R}f(x)dx}+\int_{R}^{R+i\pi}f(z)dz+\int_{R+i\pi}^{-R+i\pi}f(z)dz+\int_{-R+i\pi}^{-R}f(z)dz
$$
where the integral over $[-R,R]$ vanishes because the integrand is an odd function.
Taking $R\to\infty$ on both sides, we get
$$
\require{cancel} \lim_{R\to\infty}\oint_C f(z)dz = \cancelto{0}{\lim_{R\to\infty}\int_{R}^{R+i\pi}f(z)dz} + \lim_{R\to\infty}\int_{R+i\pi}^{-R+i\pi}f(z)dz + \cancelto{0}{\lim_{R\to\infty}\int_{-R+i\pi}^{-R}f(z)dz}\,.
$$
The integral over the vertical line segment on the right of $C$ vanishes due to the Estimation Lemma and the Squeeze Theorem, and the other integral over the left vertical segment also vanishes for the same kind of reason, I think.
To prove the integral over the vertical line segment on the right of the box decays to $0$ as $R \to \infty$, we bound the modulus of $f(z)$ first like
$$
\begin{align}
|f(z)| &= \left|\frac{z^{2}e^{2z}\left(e^{z}-1\right)^{3}\left(e^{z}+1\right)^{3}}{\left(e^{2z}+1\right)\left(6e^{4z}+e^{8z}+1\right)}\right| \\
&\leq \frac{\left|z\right|^{2}\cdot\left|e^{2z}\right|\cdot\left(3\left|e^{z}\right|+3\left|e^{2z}\right|+3\left|e^{3z}\right|+1\right)^{2}}{\left|\left|e^{2z}\right|-1\right|\cdot\left|6\left|e^{4z}\right|-\left|e^{8z}\right|-1\right|}\,. \\
\end{align}
$$
Parameterizing $z=R+iy$ for $y \in [0,\pi]$, we get
$$
\begin{align}
|f(R+iy)| &\leq \frac{\left|R+iy\right|^{2}\cdot\left|e^{2\left(R+iy\right)}\right|\cdot\left(3\left|e^{\left(R+iy\right)}\right|+3\left|e^{2\left(R+iy\right)}\right|+3\left|e^{3\left(R+iy\right)}\right|+1\right)^{2}}{\left|\left|e^{2\left(R+iy\right)}\right|-1\right|\cdot\left|6\left|e^{4\left(R+iy\right)}\right|-\left|e^{8\left(R+iy\right)}\right|-1\right|} \\
&\leq \frac{e^{2R}\left(R^{2}+2Ry+y^{2}\right)\left(3e^{R}+3e^{2R}+3e^{3R}+1\right)^{2}}{\left(e^{2R}-1\right)\left(e^{8R}-6e^{4R}+1\right)} \\
&:= M(R)\,. \\
\end{align}
$$
As $R \to \infty$, we see that $M(R)\to 0$ because the largest expression in the numerator grows slower than that of the denominator. Using the Estimation Lemma, we write
$$
0 \leq \left|\int_{R}^{R+i\pi} f(z)dz\right| \leq \pi M(R)\,.
$$
The $\pi$ comes from the length of the right side of the box contour. Taking $R \to \infty$ and employing the Squeeze Theorem, we get
$$
\lim_{R \to \infty} 0 \leq \lim_{R \to \infty} \left|\int_{R}^{R+i\pi} f(z)dz\right| \leq \pi \lim_{R \to \infty} M(R)
$$
$$
\implies \lim_{R \to \infty}\left|\int_{R}^{R+i\pi} f(z)dz\right| = 0\,.
$$
Thus,
$$
\bbox[15px,#FFF5FE]{\lim_{R \to \infty}\int_{R}^{R+i\pi} f(z)dz = 0\,.}
$$
The other integral should be evaluated in a similar fashion.
For the integral that doesn't vanish, it surprisingly helps us recover the integral we want. Parameterizing $z = x+ i\pi$, we get
$$
\begin{align}
\lim_{R\to\infty}\int_{R+i\pi}^{-R+i\pi}f(z)dz &= -\lim_{R\to\infty}\int_{-R}^{R}f(x+i\pi)d(x+i\pi) \\
&= -\lim_{R\to\infty}\int_{-R}^{R}\frac{\left(x+i\pi\right)^{2}e^{2\left(x+i\pi\right)}\left(e^{x+i\pi}-1\right)^{3}\left(e^{x+i\pi}+1\right)^{3}}{\left(e^{2\left(x+i\pi\right)}+1\right)\left(6e^{4\left(x+i\pi\right)}+e^{8\left(x+i\pi\right)}+1\right)}dx \\
\require{cancel} &= -\lim_{R\to\infty}\cancelto{0}{\int_{-R}^{R}f(x)dx} - 2\pi i \mathcal{I} + \pi^2 \lim_{R\to\infty}\cancelto{0}{\int_{-R}^{R}\frac{e^{2x}\left(e^{x}-1\right)^{3}\left(e^{x}+1\right)^{3}}{\left(e^{2x}+1\right)\left(6e^{4x}+e^{8x}+1\right)}dx}\,. \\
\end{align}
$$
The last integral vanishes — thankfully — because its integrand is an odd function. Same for the other one that canceled to $0$.
So far, we have
$$
-\frac{1}{2\pi i} \lim_{R\to\infty} \oint_C f(z)dz = \mathcal{I}\,.
$$
But we can proceed further using Cauchy's Residue Theorem and equating the real part on both sides as follows:
$$
-\Re\frac{1}{\cancel{2\pi i}} \cancel{2\pi i} \sum_{n \in \mathcal{P}}\mathop{\mathrm{Res}}_{z=n}f(z) = \Re I = I\,.
$$
Here, we will denote $\mathcal{P}$ as the set of poles inside $C$, which is $\left\{z_1(0), z_2(0), z_2(1), z_3(0), z_3(1)\right\}$.
Taking the residue at the simple pole $z_1(0)$ of $f(z)dz$, we have
$$
\begin{align}
-\Re \mathop{\mathrm{Res}}_{z=i\pi/2}f(z) &= -\Re\lim_{z \to i\pi/2}\frac{z^{2}e^{2z}\left(e^{z}-1\right)^{3}\left(e^{z}+1\right)^{3}}{\frac{d}{dz}\left(e^{2z}+1\right)\left(6e^{4z}+e^{8z}+1\right)} \\
&= -\Re\lim_{z \to i\pi/2} \frac{z^{2}e^{2z}\left(e^{z}-1\right)^{3}\left(e^{z}+1\right)^{3}}{2e^{2z}\left(12e^{2z}+18e^{4z}+4e^{6z}+5e^{8z}+1\right)} \\
&= -\Re \frac{\pi^2}{8} \\
&= -\frac{\pi^2}{8}\,. \\
\end{align}
$$
You can probably guess that the other residues at the four other elements of $\mathcal{P}$ are a drag to calculate, so I'll omit the details. Basically,
$$
\begin{align}
-\Re \mathop{\mathrm{Res}}_{z=z_2(0)}f(z) &= \frac{\pi^{2}}{128}-\frac{1}{128}\ln^{2}\left(3-2\sqrt{2}\right) \\
-\Re \mathop{\mathrm{Res}}_{z=z_2(1)}f(z) &= \frac{9\pi^{2}}{128}-\frac{1}{128}\ln^{2}\left(3-2\sqrt{2}\right) \\
-\Re \mathop{\mathrm{Res}}_{z=z_3(0)}f(z) &= \frac{\pi^{2}}{128}-\frac{1}{128}\ln^{2}\left(3+2\sqrt{2}\right) \\
-\Re \mathop{\mathrm{Res}}_{z=z_3(1)}f(z) &= \frac{9\pi^{2}}{128}-\frac{1}{128}\ln^{2}\left(3+2\sqrt{2}\right)\,, \\
\end{align}
$$
after a lot of brutal computations.
All five residues sum to
$$
\bbox[15px,border:5px groove #FF18F3 ]{-\Re\sum_{n \in \mathcal{P}}\mathop{\mathrm{Res}}_{z=n}f(z) = \frac{\pi^{2}}{32}-\frac{1}{8}\ln^{2}\left(1+\sqrt{2}\right)\,.}
$$
We finally have
$$
\bbox[15px,#FBFFEB ,border:5px groove #00906E ]{\int_{0}^{1}\frac{x^{3}\operatorname{artanh}x}{1+x^{4}}dx=\frac{\pi^{2}}{32}-\frac{1}{8}\ln^{2}\left(1+\sqrt{2}\right)}
$$
and we're done! $\blacksquare$
Cheers :)
Best Answer
Awesome observation! Rectifying your computation, let
$$ \chi_2(x) = \int_{0}^{x} \frac{\operatorname{artanh} t}{t} \, \mathrm{d}t = \sum_{k=0}^{\infty} \frac{x^{2k+1}}{(2k+1)^2} = x + \frac{x^3}{3^2} + \frac{x^5}{5^2} + \frac{x^7}{7^2} + \cdots $$
be the Legendre chi function of order $2$. Then OP's computation essentially boils down to verifying:
Once this has been established, plugging $x = 1$ gives
$$ \frac{\pi^2}{16} = \chi_2(1) - \int_{0}^{1} \frac{2t}{1+t^2} \operatorname{artanh} t \, \mathrm{d}t = \int_{0}^{1} \frac{2t}{1+t^2} \operatorname{artanh} t \, \mathrm{d}t, $$
yielding $\chi_2(1) = \frac{\pi^2}{8}$ and hence the Basel's identity $\zeta(2) = \frac{\pi^2}{6}$.
So, it remains to prove $\text{(1)}$ and $\text{(2)}$. We will closely follow OP's computation (but in a rigorous manner) in order to establish both identites.
Proof of $(1)$. For $|x| < 1$, expanding the integrand as a power series in $t$ and applying Fubini's theorem gives
\begin{align*} -\int_{0}^{x} \frac{2t}{1+t^2} \operatorname{artanh}(xt) \, \mathrm{d}t &= -2 \int_{0}^{x} \biggl( \sum_{j=0}^{\infty} (-1)^j t^{2j+1} \biggr) \biggl( \sum_{k=0}^{\infty} \frac{(xt)^{2k+1}}{2k+1} \biggr) \, \mathrm{d}t \\ &= -2 \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} (-1)^j \frac{x^{2k+1}}{2k+1} \int_{0}^{x} t^{2j+2k+2} \, \mathrm{d}t \\ &= 2 \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \frac{(-1)^k x^{2k+1}}{2k+1} \cdot \frac{(-1)^{k+j+1}x^{2k+2j+3}}{2k+2j+3} \\ &= 2 \sum_{0\leq k < l} \frac{(-1)^k x^{2k+1}}{2k+1} \cdot \frac{(-1)^{l}x^{2l+1}}{2l+1} \\ &= \biggl( \sum_{k=0}^{\infty} \frac{(-1)^k x^{2k+1}}{2k+1} \biggr)^2 - \sum_{k=0}^{\infty} \frac{x^{2(2k+1)}}{(2k+1)^2} \\ &= (\arctan x)^2 - \chi_2(x^2). \end{align*}
This proves $\text{(1)}$ for $|x| < 1$. Since both sides of $\text{(1)}$ are continuous at $x = 1$ from the left, the validity of $\text{(1)}$ extends to $x = 1$, proving the desired identity.
Proof of $(2)$. Let $|x| < 1$. Proceeding similarly as above,
\begin{align*} \int_{0}^{x} \frac{2t}{1+t^2} \operatorname{artanh}(t/x) \, \mathrm{d}t &= \int_{0}^{1} \frac{2x^2t}{1+x^2t^2} \operatorname{artanh} t \, \mathrm{d}t \\ &= 2 \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \frac{(-1)^j}{2k+1} x^{2j+2} \int_{0}^{1} t^{2k+2j+2} \, \mathrm{d}t \\ &= 2 \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \frac{(-1)^j x^{2j+2}}{(2k+1)(2k+2j+3)} \\ &= \sum_{j=0}^{\infty} (-1)^j x^{2j+2} \biggl( \sum_{k=0}^{\infty} \frac{2}{(2k+1)(2k+2j+3)} \biggr) \end{align*}
Now we focus on the inner sum in the last line. Applying partial fraction decomposition,
\begin{align*} \sum_{k=0}^{\infty} \frac{2}{(2k+1)(2k+2j+3)} &= \frac{2}{2j+2} \sum_{k=0}^{\infty} \biggl( \frac{1}{2k+1} - \frac{1}{2k+2j+3} \biggr) \\ &= \frac{2}{2j+2} \sum_{k=0}^{j} \frac{1}{2k+1} \\ &= \frac{1}{2j+2} \biggl(\sum_{k=0}^{j} \frac{1}{2k+1} + \sum_{k=0}^{j} \frac{1}{2(j-k)+1} \biggr) \\ &= \sum_{k=0}^{j} \frac{1}{(2k+1)(2(j-k)+1)}. \end{align*}
Plugging this back,
\begin{align*} \int_{0}^{x} \frac{2t}{1+t^2} \operatorname{artanh}(t/x) \, \mathrm{d}t &= \sum_{j=0}^{\infty} \sum_{k=0}^{j} \frac{(-1)^k x^{2k+1}}{2k+1} \cdot \frac{(-1)^{j-k} x^{2(j-k)+1}}{2(j-k)+1}, \end{align*}
which is the Cauchy product of the power series of $\arctan x$ with itself, completing the proof.