We can think of $I_{n}$ as being a classical partition function for $n$ beads on a circle which cannot pass through each other, with logarithmic interaction potential between each bead and its next-to-nearest neighbors on either side. For $I_{2n}$ the beads fall into two ``colors" which do not have logarithmic interactions with each other; while for $I_{2n+1}$ the beads do not fall into two independent groups.
We make two changes of variable. First, we can label the coordinates of the $k^{th}$ bead as $y_k$, where $y_1=0$ is fixed (exploiting the translation invariance of the problem) and we define $y_{2n+k} = 1+y_k$ (because of the periodic nature of the circle):
$$
u_i = y_{i+1}-y_i\ ,\qquad y_1=0\ ,\qquad y_{2n+i}\equiv 1+y_i \ .
$$
Then the integral can be written as a path ordered expression without the delta function constraint as
$$
I_{n}= \int_0^1 dy_{n} \int_0^{y_{n}} dy_{n-1}\cdots\int_0^{y_3} dy_2\, \prod_{k=1}^{2n}\frac{1}{y_{k+2}-y_k}\ .
$$
The second change of variables to $\{y_2,\ldots,y_n\}\to \{s_2,\ldots,s_n\}$ in order change the integration domain to a unit hypercube:
$$
y_{k} =\prod_{j=k}^{n} s_{j}\ ,
$$
with Jacobian
$$
J_n = \prod_{j=3}^{n} s_j^{j-2}\ .
$$
With this change of variables, $I_{2n}$ becomes (for $n\ge 2$)
$$
I_{2n} = \int_0^1 d^{2n-1}{\bf s}\, \prod_{j=2}^{2n-1} \,
\frac{1}{1-s_j s_{j+1}} \frac{1}{1-s_{2n}{\bf S}_{2n+1}}\equiv \int_0^1d^{2n-1}{\bf s}\, {\cal F}_{2n}({\bf s})\ $$
where $d^{2n-1}{\bf s}=ds_2\cdots ds_{2n}$, and the integral sign indicates that each of the $s$ variables is being integrated from zero to one, and we have defined
$$
{\bf S}_{k}\equiv 1+(-1)^k\prod_{j=2}^{k-2} s_j\ ,\qquad
{\cal F}_{2n}({\bf s}) = \prod_{j=2}^{2n-1} \,
\frac{1}{1-s_j s_{j+1}} \frac{1}{1-s_{2n}{\bf S}_{2n+1}}\ ,
$$
with
$$
S_3=0\ ,\qquad {\cal F}_{2}({\bf s}) =1\ .
$$
Note that for odd $k$, ${\bf S}_k<1$, while for even $k$, ${\bf S}_k>1$.
This object ${\bf S}_k$ has the property for any $k$
$$
{\bf S}_{k+1} -s_{k-1} = 1-s_{k-1} {\bf S}_k\ .
$$
The strategy is to consider developing a recursion relation when integrating over $ds_{2n}$ and $ds_{2n-1}$, relating $I_{2n}$ to $I_{2n-2}$. To that end it is useful to define the following functions of $x$, $y$ in the domain $ 0<x<1,\ 0<y<1$:
$$
{\cal P}_k(x,y) = \frac{1}{(2k)!}
\prod_{i=1}^k
\left(\pi^2 (2k-1)^2 + \ln^2\left[\frac{1-x}{x(1-y)}\right]\right)\ ,\qquad {\cal P}_0(x,y)\equiv 1\ ,
$$
and
$$
{\cal G}(\alpha,x,y) = \sum_{n=0}^\infty \left(\frac{\alpha}{\pi}\right)^{2n} {\cal P}_n(x,y) =
\frac{1}{2 \sqrt{1-\alpha ^2}}\left[\left(\frac{1-x}{x(1-y)}\right)^{c}+\left(\frac{1-x}{x(1-y) }\right)^{-c}\,\right]\ ,
$$
$$
c\equiv \frac{\sin ^{-1}(\alpha )}{\pi }\ .
$$
We generalize the problem to considering the integral
$$
{\cal I}_{2n}(\alpha) = \int_0^1d^{2n-1}{\bf s}\, {\cal F}_{2n}({\bf s})\,{\cal G}(\alpha,s_{2n},{\bf S}_{2n+1})\ .
$$
We can perform the $s_{2n}$ and $s_{2n-1}$ integrals in ${\cal I}_{2n}$ using the results (using the properties of ${\bf S_k}$ above)
For $0<s_{2n-1}<1$ and $0<{\bf S}_{2n+1}<1$:
$$
\frac{1}{2} \int_0^1 ds_{2n} \frac{1}{(1-s_{2n-1} s_{2n})(1-s_{2n}{\bf S}_{2n+1})}\left[ \left(\frac{1-s_{2n}}{s_{2n}(1-{\bf S}_{2n+1})}\right)^c+ \left(\frac{1-s_{2n}}{s_{2n}(1-{\bf S}_{2n+1})}\right)^{-c}\right]
= \frac{\pi \csc (\pi c) \left(\left(\frac{1-s_{2n-1}}{1-{\bf S}_{2n+1}}\right)^{-c}-\left(\frac{1-s_{2n-1}}{1-{\bf S}_{2n+1}}\right)^c\right)}{2(s_{2n-1}-{\bf S}_{2n+1})}
=\frac{\pi \csc (\pi c)\left(\left(\frac{1-s_{2n-1}}{s_{2n-1}({\bf S}_{2n}-1)}\right)^{-c}-\left(\frac{1-s_{2n-1}}{s_{2n-1}({\bf S}_{2n}-1)}\right)^c\right)}{2(1-s_{2n-1}{\bf S}_{2n})}\
$$
For $0<s_{2n-2}<1$ and $1<{\bf S}_{2n}$:
$$
\frac{1}{2} \int_0^1 ds_{2n-1} \frac{1}{(1-s_{2n-2}s_{2n-1})(1-{\bf S}_{2n}s_{2n-1})}
\left[ \left(\frac{1-s_{2n-1}}{s_{2n-1}({\bf S}_{2n}-1)}\right)^c- \left(\frac{1-s_{2n-1}}{s_{2n-1}({\bf S}_{2n}-1)}\right)^{-c}\right]
= -\frac{\pi \csc (\pi c)
\left(\left(\frac{1-s_{2n-2}}{{\bf S}_{2n}-1}\right)^{-c}+\left(\frac{1-s_{2n-2}}{{\bf S}_{2n}-1}\right)^c-2 \cos
(\pi c)\right)}{2 (s_{2n-2}-{\bf S}_{2n})}
= -\frac{\pi \csc (\pi c)
\left(\left(\frac{1-s_{2n-2}}{1-s_{2n-2}{\bf S}_{2n-1}}\right)^{-c}+\left(\frac{1-s_{2n-2}}{1-s_{2n-2}{\bf S}_{2n-1}}\right)^c-2 \cos
(\pi c)\right)}{2 (1-s_{2n-2}{\bf S}_{2n-1})}
$$
With these integrals we can perform the integrations over $s_{2n}$ and $s_{2n-1}$ in our generalized integral ${\cal I}_{2n}(\alpha)$, obtaining
$$
{\cal I}_{2n}(\alpha)
=\int d^{2n-3}{\bf s} \, {\cal F}_{2n-2}({\bf s})\, \left[\pi^2\csc^2(c\pi)\left({\cal G}(\alpha,s_{2n-2},{\bf S}_{2n-1}) -\frac{ \cos c\pi}{ \sqrt{1-\alpha^2}}\right) \right]
\,
=\int d^{2n-3}{\bf s}{\cal F}_{2n-2}({\bf s})\, \left[\frac{\pi^2}{\alpha^2}\left({\cal G}(\alpha,s_{2n-2},{\bf S}_{2n-1})-1\right) \right]
$$
where to get the second line we just plugged in $\pi c=\sin^{-1}\alpha$. Referring to the definition of ${\cal G}$ in eq.(\ref{gdef}), we can equate powers of $\alpha$ on both sides of the above equation with the result that for every $k\ge 0$,
$$
\int_0^1d^{2n-1}{\bf s}\, {\cal F}_{2n}({\bf s})\, {\cal P}_k(s_{2n},{\bf S}_{2n+1})
= \int_0^1d^{2n-3}{\bf s}\, {\cal F}_{2n-2}({\bf s})\, {\cal P}_{k+1}(s_{2n-2},{\bf S}_{2n-1})
$$
which is a pretty result.
The above result allows us to write for the desired $2n$-dimensional integrals as one-dimensional integrals
$$
I_{2n}=\int_0^1d^{2n-1}{\bf s}\, {\cal F}_{2n}({\bf s})\, {\cal P}_0(s_{2n},{\bf S}_{2n+1})
=\int_0^1ds_2 {\cal P}_{n-1}(s_{2},0)\ .
$$
The above results then imply that
$$
\sum_{n=0}^\infty \left(\frac{\alpha}{\pi}\right)^{2n} I_{2n+2} = \int_0^1dx\,\sum_{n=0}^\infty \left(\frac{\alpha}{\pi}\right)^{2n} {\cal P}_n(x,0)
= \int_0^1dx\, {\cal G}(\alpha,x,0)
=\int_0^1dx\frac{1}{2 \sqrt{1-\alpha ^2}}\left[\left(\frac{1-x}{x}\right)^{c}+\left(\frac{1-x}{x }\right)^{-c}\right]
= \frac{\sin^{-1}\alpha}{\alpha\sqrt{1-\alpha^2}}
= \sum_{n=0}^\infty (2\alpha)^{2n} B(n+1,n+1)\ .
$$
Equating powers of $\alpha$ between the first and last expressions answers the posted question.
This solution was found in collaboration with E. Mereghetti (we're physicists, so the language might look odd).
