Just some remarks from a higher perspective. Integrals:
$$\tag{1} \int_0^1 f(x) dx$$
when $f(x)$ is one of
$$\small \frac{\log^n x \log(1+x^2)}{1\pm x},\frac{\log^{2n} x \log(1+x^2)}{1+x^2},\frac{x\log^{2n} x \log(1+x^2)}{1+x^2},\frac{\log^{2n} x \arctan x}{1- x},\frac{\log^{n} x \arctan x}{1+ x},\frac{x\log^{2n} x \arctan x}{1 + x^2},\frac{\log^{2n+1} x \arctan x}{1 + x^2},\frac{\log^{2n} x\log(1\pm x)}{1 + x^2},\frac{x\log^{2n+1}x \log(1\pm x)}{1 + x^2},\frac{\log^{2n}x \log(1\pm x) \log(1+x^2)}{x},\frac{\log^{2n+1}x \arctan x \log(1+x^2)}{x},\frac{\log^{2n+1}x \arctan x \log(1\pm x)}{x},\frac{\log^{2n+1} x \text{Li}_{2m}(x)}{1+x^2},
\frac{\log^{2n} x \text{Li}_{2m+1}(x)}{1+x^2},
\frac{x \log^{2n+1} x \text{Li}_{2m+1}(x)}{1+x^2} \text{ etc...}$$
all have results of form
$$\tag{2}\int_0^1 f(x) dx = \sum_{i=1}^{w} a_i L(i)L(w-i)$$
here $a_i \in \mathbb{Q}$, $L(s) =\zeta(s)$ or $\beta(s)$ are Riemann zeta or Dirichlet beta, $w$ is the weight of integrand, namely degree of numerator plus one.
OP's $I_2$ belongs to above cases; $I_1$ also, had the power of $\log$ were even; the odd power case is more difficult.
For $I_1$, the fact that MZIntegrate can do it already implies a mechanical proof exists. If you want magical/elegant solutions, perhaps Conrad Valean or Ali Shadhar can shed light on this. A highly non-trivial extension is $\int_{0}^{1}\frac{\log^3{(x)}\log{(1\pm x)}\log{(1+x^2)}}{x}dx$, you can put it into MZIntegrate to see what it looks like, I think Conrad and Ali have yet to found a solution to this one.
For example:
$$\int_0^1 \frac{\log ^5(x) \log \left(x^2+1\right) \arctan x}{x} dx = 720 \beta(8)-\frac{25 \pi ^5 \zeta (3)}{32}-\frac{15 \pi ^3 \zeta (5)}{2}-\frac{123825 \pi \zeta (7)}{2048}$$
$$\int_0^1 \frac{\log^4 x \log(1+x)}{1+x^2} dx = \frac{7 \pi ^4 G}{30}+2 \pi ^2 \beta(4)-48 \beta(6)+\frac{5}{128} \pi ^5 \log (2)$$
$$\int_0^1 \frac{x \log^5 x \log(1-x)}{1+x^2} dx = -\frac{\pi ^4 \zeta (3)}{8}-\frac{75 \pi ^2 \zeta (5)}{128}+\frac{129495 \zeta (7)}{2048}-\frac{33}{512}\pi ^6 \log (2)$$
$$\int_0^1 \frac{\log^4 x \log(1-x) \log(1+x^2)}{x} dx = -\frac{5 \pi ^5 G}{32}-\frac{1}{2} 3 \pi ^3 \beta(4)-12 \pi \beta(6)+\frac{121 \pi ^4 \zeta (3)}{2560}+\frac{105 \pi ^2 \zeta (5)}{512}+\frac{121557 \zeta (7)}{1024}$$
when weight $w\geq 7$, MZIntegrate is not able to evaluate them directly. I will add this functionality on the next update.
If one replaces $1+x^2$ occurring in $f(x)$ by $1+x+x^2$ or $1-x+x^2$ or $1+x^4$ etc, one obtain similar expressions, except now $L$ has to be other Dirichlet L-functions.
Proof sketch: for all cases $f(x)$ above, the integral $\int_0^1 f(x) dx$ can be converted to colored multiple zeta values (CMZVs) of depth 2 and level 4 and weight $w$. There is a result saying that whenever depth and weight of a CMZV satisfies certain parity condition, its real/imaginary part reduces to lower depths CMZVs. In our case, we have reduction into depth 1 CMZVs, which when level = 4 is simply $\zeta(s)$ or $\beta(s)$, giving $(2)$. QED.
MZIntegrate
is just a command within that package. – heropup Jul 19 '22 at 09:22