With the help of programs I have been able to conjecture $$\int _0^1\int _0^z\int _0^y\frac{1}{\left(1-x^2\right)\left(1+y^2\right)\left(1+z^2\right)}\:\mathrm{d}x\:\mathrm{d}y\:\mathrm{d}z=\frac{\pi }{8}G-\frac{7}{32}\zeta \left(3\right)$$ But is there a simple way to prove this?
What I have got thus far is: $$\int _0^1\int _0^z\int _0^y\frac{1}{\left(1-x^2\right)\left(1+y^2\right)\left(1+z^2\right)}\:\mathrm{d}x\:\mathrm{d}y\:\mathrm{d}z=-\frac{1}{2}\int _0^1\frac{1}{1+z^2}\int _0^z\frac{\ln \left(\frac{1-y}{1+y}\right)}{1+y^2}\:\mathrm{d}y\:\mathrm{d}z$$ $$=-\frac{1}{2}\int _0^1\frac{\ln \left(\frac{1-y}{1+y}\right)}{1+y^2}\int _y^1\frac{1}{1+z^2}\:\mathrm{d}z\:\mathrm{d}y=-\frac{1}{2}\int _0^1\frac{\left(\frac{\pi }{4}-\arctan \left(y\right)\right)\ln \left(\frac{1-y}{1+y}\right)}{1+y^2}\:\mathrm{d}y$$ But I'm not sure how to advance from here, is the current path I'm taking correct?