How can we reproduce the exact solution of the integral of the splitting function, which play a role in the dipole factorisation formulae?

Question

I have a problem performing the following integration provided in the paper by Catani and Seymour (https://arxiv.org/abs/hep-ph/9605323) page 27.

Given is the integral

$$
\mathcal{V}=\int_0^1 (z(1-z))^{-\epsilon} \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}V(z;y) dydz
$$

with

$$
V(z;y)=\frac{2}{1-z(1-y)}-(1+z)-\epsilon(1-z).
$$

In the paper the exact result and an approximation is stated as follows:

$$
\mathcal{V}=\frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon^2}+\frac{1}{\epsilon}\frac{3+\epsilon}{2(1-3\epsilon)}\right]=\frac{1}{\epsilon^2}+\frac{3}{2\epsilon}+5-\frac{\pi^2}{2}+\mathcal{O}(\epsilon).
$$

To me it is obvious how to deal with the second and the third summand of $V$ under the integral. I simply expand and employ the definition of the Euler Beta function

$$
B(a,b)=\int_0^1 t^{a-1}(1-t)^{b-1} dt=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}.
$$

Hence

\begin{align}
\int_0^1 (z(1-z))^{-\epsilon}& \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}\left[-(1+z)-\epsilon(1-z)\right] dydz \\
& = \frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon}\frac{3+\epsilon}{2(1-3\epsilon)}\right]
\\ & = \frac{3}{2\epsilon}+5+\mathcal{O}(\epsilon)
\end{align}

for the simple part.
Do you have any suggestions how to treat the singular term of $V$ in order to get the exact result? In particular how to solve

\begin{align}
\int_0^1 (z(1-z))^{-\epsilon}& \int_0^1 (1-y)^{1-2\epsilon}y^{-1-\epsilon}\left[\frac{2}{1-z(1-y)}\right] dydz\\
&=\frac{\Gamma(1-\epsilon)^3}{\Gamma(1-3\epsilon)}\left[\frac{1}{\epsilon^2}\right]\\
&=\frac{1}{\epsilon^2}-\frac{\pi^2}{2}+\mathcal{O}(\epsilon).
\end{align}

Thank you in advance.

Edit: After some effort in reverse engineering (basically rewriting the product of Gamma-functions in the exact result in terms of two Euler Beta functions) the initial problem boils down to showing the equality of

\begin{align}
& \int_0^1 (z(1-z))^{-\epsilon} \int_0^1 y^{1-2\epsilon}(1-y)^{-1-\epsilon}\left[\frac{2}{1-zy}\right] dydz \\
&=\int_0^1 (t(1-t))^{-\epsilon} \int_0^1 s^{1-2\epsilon}(1-s)^{-1-\epsilon}\left[\frac{(1-s)}{t(1-t)s^2}\right] dsdt.
\end{align}

I guess the solution is an integral transformation. Does anyone know how to transform in order to show the equality of the upper expressions?