Cross-posting this integral on AoPS brings Y. Sharifi's solution here after a day. Quite amazing one!
I will copy here his entire solution:
Let $I$ be your integral. Using the identity $\ln x \ln(1-x)+\text{Li}_2(x)=\zeta(2)-\text{Li}_2(1-x),$ we have
$$I=-\int_0^1 \frac{\text{Li}_2(x)\text{Li}_2(1-x)}{x(1-x)} \ dx + \zeta(2)\int_0^1 \left(\frac{\text{Li}_2(x)}{x(1-x)}-\frac{\zeta(2)}{1-x}+\frac{\text{Li}_2(1-x)}{1-x}\right)dx.$$
Let
$$J=\int_0^1 \frac{\text{Li}_2(x)\text{Li}_2(1-x)}{x(1-x)} \ dx, \ \ \ \ \ K:=\int_0^1 \left(\frac{\text{Li}_2(x)}{x(1-x)}-\frac{\zeta(2)}{1-x}+\frac{\text{Li}_2(1-x)}{1-x}\right)dx.$$
So
$$I=\zeta(2)K - J. \ \ \ \ \ \ \ \ \ (1)$$
We first show that $K=0.$ Start with using integration by parts in $K,$ with $u=\frac{\text{Li}_2(x)}{x}-\zeta(2)+\text{Li}_2(1-x)$ and $dv=\frac{dx}{1-x}.$ Then
$$K=\int_0^1 \ln(1-x)\left(\frac{\ln x}{1-x}-\frac{\ln(1-x)}{x^2}-\frac{\text{Li}_2(x)}{x^2}\right)dx. \ \ \ \ \ \ \ \ \ \ (2)$$
Using the Maclaurin series of $\ln(1-x),$ we quickly find the first integral in $K$
$$\int_0^1 \frac{\ln x \ln(1-x)}{1-x} \ dx = \int_0^1 \frac{\ln x \ln(1-x)}{x} \ dx=\zeta(3). \ \ \ \ \ \ \ \ \ \ (3)$$
Next, we ignore the second integral in $K$ for now and we look at the third one, i.e. $\int_0^1 \frac{\ln(1-x) \text{Li}_2(x)}{x^2} \ dx.$ In this integral, we use integration by parts with $u=\ln(1-x)\text{Li}_2(x)$ and $dv=\frac{dx}{x^2};$ notice that we need to choose $v=1-\frac{1}{x}.$ So
$$\int_0^1 \frac{\ln(1-x) \text{Li}_2(x)}{x^2} \ dx=\int_0^1\left(1-\frac{1}{x}\right)\left(\frac{\text{Li}_2(x)}{1-x}+\frac{\ln^2(1-x)}{x}\right) dx$$
$$=-\int_0^1 \frac{\text{Li}_2(x)}{x} \ dx + \int_0^1 \frac{\ln^2(1-x)}{x} \ dx - \int_0^1 \frac{\ln^2(1-x)}{x^2} \ dx=-\zeta(3)+\int_0^1 \frac{\ln^2x}{1-x} \ dx -\int_0^1 \frac{\ln^2(1-x)}{x^2} \ dx.$$
$$=\zeta(3)-\int_0^1 \frac{\ln^2(1-x)}{x^2} \ dx. \ \ \ \ \ \ \ \ \ (4)$$
Thus, by $(2),(3)$ and $(4),$ we have $K=0$ and hence, by $(1),$
$$I=-J=-\int_0^1 \frac{\text{Li}_2(x)\text{Li}_2(1-x)}{x(1-x)} \ dx=-2\int_0^1 \frac{\text{Li}_2(x)\text{Li}_2(1-x)}{x} \ dx.$$
So integration by parts with $u=\text{Li}_2(1-x)$ and $dv=\frac{\text{Li}_2(x)}{x} \ dx$ gives
$$I=2\int_0^1 \frac{\text{Li}_3(x)\ln x}{1-x} \ dx=2\int_0^1 \text{Li}_3(x) \ln x \sum_{m \ge 1}x^{m-1} dx=2\sum_{m \ge 1} \int_0^1 x^{m-1}\text{Li}_3(x) \ln x \ dx$$
$$=2\sum_{m \ge 1} \int_0^1x^{m-1}\sum_{n \ge 1} \frac{x^n}{n^3} \ln x \ dx=2\sum_{m,n \ge 1} \frac{1}{n^3}\int_0^1x^{n+m-1}\ln x \ dx=-2\sum_{m,n \ge 1} \frac{1}{n^3(n+m)^2}$$
$$=-\sum_{m,n \ge 1} \left(\frac{1}{n^3(n+m)^2}+\frac{1}{m^3(n+m)^2}\right). \ \ \ \ \ \ \ \ \ (5)$$
So $(5)$ and the following identity
$$\frac{1}{n^3(n+m)^2}+\frac{1}{m^3(n+m)^2}=\frac{1}{n^3m^2}-\frac{2}{n^2m^3}+\frac{3}{m^3n(n+m)}$$
together give
$$I=-\sum_{m,n \ge 1}\left(\frac{1}{n^3m^2}-\frac{2}{n^2m^3}+\frac{3}{m^3n(n+m)}\right)=\zeta(2)\zeta(3)-3\sum_{m,n \ge 1} \frac{1}{m^3n(n+m)}$$
$$=\zeta(2)\zeta(3)-3\sum_{m \ge 1} \frac{1}{m^4} \sum_{n \ge 1}\left(\frac{1}{n}-\frac{1}{n+m}\right)=\zeta(2)\zeta(3)-3\sum_{m \ge 1} \frac{H_m}{m^4}, \ \ \ \ \ \ \ \ \ (6)$$
where, as usual, $H_m:=\sum_{j=1}^m \frac{1}{j}$ is the $m$-th harmonic number. Now we use Euler's formula
$$\sum_{m \ge 1} \frac{H_m}{m^k}=\left(1+\frac{k}{2}\right)\zeta(k+1)-\frac{1}{2}\sum_{i=1}^{k-2}\zeta(i+1)\zeta(k-i), \ \ \ \ k \ge 2,$$
with $k=4$ to get
$$\sum_{m \ge 1} \frac{H_m}{m^4}=3\zeta(5)-\zeta(2)\zeta(3)$$
and so, by $(6),$
$$I=\zeta(2)\zeta(3)-3(3\zeta(5)-\zeta(2)\zeta(3))=4\zeta(2)\zeta(3)-9\zeta(5).$$
Edit.
This integral was proposed two years ago in RMM and it appeared as problem UP $089$.
See in this link, at the page $70$.
Assuming that $s$ is an integer $\geq 2$, each term in the binomial expansion of the integrand is separable into a product of integrals over $x$ and $y$. Defining $N=s-2$ and $z'=1-z$ for notational simplicity,
$$
I(1-z',N+2) = \int_{0}^1 dx \int_{0}^1 dy \sum_{k=0}^N \binom{N}{k} \left( -\frac{(1-x)(1-y)}{(1-z'x)(1-z'y)}\right)^k \\
= \sum_{k=0}^N \binom{N}{k} (-1)^k \left\{
\left(
\int_{0}^1 \left( \frac{1-x}{1-z'x}\right)^k \, dx
\right)^2
\right\}
$$
As a check, when $z=1$ and thus $z'=0$, the inner integral is $1/(1+k)$, and so
$$
I(1,N+2) = \sum_{k=0}^N \binom{N}{k} \frac{(-1)^k}{(1+k)^2}
$$
which sums (as verified by Mathematica) to the same result derived by you. For general $z'$, define
$$
\xi(z) =\int_0^1 \left( \frac{1-x}{1-(1-z)x} \right)^k \, dx
$$
Unfortunately I cannot evaluate this integral, but it can be expanded to linear order around $z=1$ to give
$$
\xi(z) = \frac{1}{1+k} + \frac{k}{(k+1)(k+2)} (1-z) + \mathcal{O}((1-z)^2)
$$
yielding (for integral $s \geq 2$ at least) (edit and with some corrections to converting the expansion of $\xi(z)$ into that of $\xi^2(z)$)
$$
I(z,s) = \frac{\psi(s)+\gamma}{s-1} + 2 \frac{s(\psi(s+1)+\gamma)-2s+1}{s(s-1)} (1-z) + \mathcal{O}((1-z)^2)
$$
ADDITION
As pointed out by the OP, from Mathematica or integral tables (c.f. https://dlmf.nist.gov/15.6#E1) we have
$$
\xi(z) = \int_0^1 \left(\frac{1-x}{1-(1-z)x}\right)^k \, dx = \frac{1}{k+1} F(1,k;k+2;1-z)
$$
(Note that the referenced table represents the identity using the Olver form of the hypergeometric function $\mathbf{F}(a,b;c;z)=F(a,b;c;z)/\Gamma(c)$, and that $F(a,b; c; z)$ is symmetric in $a$ and $b$). This gives the result, for integral $s \geq 2$,
$$
\left. I(z,s)\right|_{s \in \mathbb{N}, s\geq 2} = \sum_{k=0}^{s-2} \binom{s-2}{k} \frac{(-1)^k}{(1+k)^2} F(1,k;k+2,1-z)^2
$$
At least formally, without worrying about convergence, this could be generalized to non-integer $s$ (and $s=0,1$) by extending the series to infinity (with generalized binomial coefficients), but this is probably less useful. The Taylor series of the squared hypergeometric function is
$$
F(1,k;k+2;1-z)^2 = \left( \sum_{n=0}^\infty \frac{(1)_n (k)_n}{(k+2)_n} \frac{(1-z)^n}{n!} \right)^2 \\
= \left( \sum_{n=0}^\infty \frac{k(k+1)}{(k+n)(k+n+1)} (1-z)^n \right)^2 \\
= \sum_{n=0}^\infty \left\{ \sum_{m=0}^n \frac{k(k+1)}{(k+m)(k+m+1)} \frac{k(k+1)}{(k+n-m)(k+n-m+1)} \right\} (1-z)^n
$$
This inner sum is somewhat messy, but can be broken up with partial fractions to yield
$$
\sum_{n=0}^\infty \frac{2k^2(k+1)^2}{2k+n+1} \left( \frac{H_{k+n}-H_{k-1}}{2k+n} - \frac{H_{k+n+1}-H_k}{2k+n+2} \right) (1-z)^n
$$
Substituting into the expression for $I(s,z)$ and exchanging the order of summation, the Taylor series $I(z,s) = \sum_{n=0}^\infty C_n (1-z)^n$ has coefficients
$$
C_n = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+n+1} \left( \frac{H_{k+n}-H_{k-1}}{2k+n} - \frac{H_{k+n+1}-H_k}{2k+n+2} \right)
$$
This reduces to the special cases above for $n=0,1$, e.g.
$$
C_0 = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+1} \left( \frac{1}{2k^2} - \frac{1}{2(k+1)^2} \right) = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{k^2 (-1)^k}{(k+1)^2} \\
C_1 = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k^2 (-1)^k}{2k+2} \left( \frac{1}{k(k+1)} - \frac{1}{(k+1)(k+2)} \right) = \sum_{k=0}^{\infty} \binom{s-2}{k} \frac{2k (-1)^k}{(k+1)^2(k+2)}
$$
and so on.
Best Answer
A similar problem and solution can be found here. Proposed by Khalef and solved by Sujee.
Since $\int_0^1 \frac{dx}{1+x^2}=\frac{\pi}{4}$, we have
$$\frac{\pi^2}{16}=\int_0^1\int_0^1\frac{dydx}{(1+x^2)(1+y^2)}\overset{t=xy}{=}\int_0^1\int_0^x\frac{dtdx}{x(1+x^2)(1+t^2/x^2)}$$
$$=\frac12\int_0^1\int_t^1\frac{dxdt}{x(1+x^2)(1+t^2/x^2)}\overset{x^2\to x}{=}\frac12\int_0^1\left(\int_{t^2}^1\frac{dx}{(1+x)(x+t^2)}\right)dt$$
$$=-\frac12\int_0^1\frac{\ln\left(\frac{4t^2}{(1+t^2)^2}\right)}{1-t^2}dt\overset{t=\frac{1-x}{1+x}}{=}-\frac12\int_0^1\frac{\ln\left(\frac{1-x^2}{1+x^2}\right)}{x}dx$$
$$\overset{x^2\to x}{=}-\frac14\int_0^1\frac{\ln\left(\frac{1-x}{1+x}\right)}{x}dx=-\frac14\int_0^1\frac{\ln\left(\frac{(1-x)^2}{1-x^2}\right)}{x}dx$$
$$=-\frac12\int_0^1\frac{\ln(1-x)}{x}dx+\frac14\underbrace{\int_0^1\frac{\ln(1-x^2)}{x}dx}_{x^2\to x}$$
$$=-\frac38\int_0^1\frac{\ln(1-x)}{x}dx\Longrightarrow \int_0^1\frac{-\ln(1-x)}{x}dx=\frac{\pi^2}{6}$$
Remark:
This solution can be considered a proof that $\zeta(2)=\frac{\pi^2}{6}$ as we have $\int_0^1\frac{-\ln(1-x)}{x}dx=\text{Li}_2(x)|_0^1=\text{Li}_2(1)=\zeta(2)$