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).
Part of what makes this question subtle is that what's intuitive depends on your background knowledge. In particular, the question of what counts as an "explicit tiny quantity" is hard to pin down. For example, Bailey, Borwein, and Borwein pointed out the example
$$
\left(\frac{1}{10^5} \sum_{n=-\infty}^\infty e^{-n^2/10^{10}}\right)^2 \approx \pi,
$$
which holds to over 42 billion digits. If you know something about Poisson summation, it's pretty obvious that this is a fraud: Poisson summation converts this series into an incredibly rapidly converging series with first term $\pi$, after which everything else is tiny. Even if you don't know about Poisson summation, the $10^5$ and $10^{10}$ should raise some suspicions. However, it's a pretty amazing example if you haven't seen this sort of thing before.
As for the third question, the truth is more dramatic than that. There are finite identities that are independent of ZFC, or any reasonable axiom system, so you don't even need anything fancy like infinite sums or integrals. This follows from a theorem of Richardson (D. Richardson, Some Undecidable Problems Involving Elementary Functions of a Real Variable, Journal of Symbolic Logic 33 (1968), 514-520, http://www.jstor.org/stable/2271358).
Specifically, consider the expressions obtainable by addition, subtraction, multiplication, and composition from the initial functions $\log 2$, $\pi$, $e^x$, $\sin x$, and $|x|$. Richardson proved that there is no algorithm to decide whether such an expression defines the zero function. Now, if the function is not identically zero then it can be proved not to be zero (find a point at which it doesn't vanish and compute it numerically with enough accuracy to verify that it is nonzero). If all the identically zero cases could be proved in ZFC too, then that would give a decision procedure: do a brute force search through all possible proofs until you find one that resolves the question. Thus, there exists an expression that actually is identically zero but where there is no proof in ZFC. In fact, if you go through Richardson's proof you can explicitly construct such an expression, although it will be terrifically complicated. (To be precise, you'll be able to prove in ZFC that if ZFC is consistent, then this identity is unprovable, but by Gödel's second incompleteness theorem you won't be able to prove that ZFC is consistent in ZFC.)
Being independent of ZFC doesn't mean these identities are true for no reason, and in fact Richardson's proof gives a reason. However, his paper shows that there is no systematic way to test identities. You can indeed end up stuck in a situation where a proposed identity really looks true numerically but you just can't find a proof.
Best Answer
The answer is 'no'. Making the substitution $$ x = \frac{(t-1)(t-5)(t^2+2t+5)}{16t^2}, $$ one finds $$ {\textstyle\sqrt{x+\sqrt{x+\sqrt{x+1}}}\,\mathrm{d}x} = \frac{(t^2-2t+5)(t^2-5)\sqrt{t^4{-}2t^2{-}40t+25}\ \mathrm{d}t}{32t^4}. $$ Denote the right hand side of the above equation by $\beta$. Now, setting $$ Q(t) = \frac{(25-25t-14t^2-3t^3+t^4)}{96t^3}, $$ one finds that $$ \beta = \mathrm{d}\left(Q(t)\sqrt{t^4{-}2t^2{-}40t+25}\right) + \frac{(10+15t-2t^2+t^3)\,\mathrm{d}t}{8t \sqrt{t^4-2t^2-40t+25}}\, $$ The second term on the right hand side is in normal form for differentials on the elliptic curve $s^2 = t^4-2t^2-40t+25$ that are odd with respect to the involution $(t,s)\mapsto(t,-s)$. Since this term represents a differential that has two poles of order $2$ (over the two points where $t=\infty$) and two poles of order $1$ (over the two points where $t=0$), an application of Liouville's Theorem (on integration in elementary terms, with the differential field taken to be the field of meromorphic functions on the elliptic curve) shows that this term is not integrable in elementary terms. Hence $\beta$ is not, and hence the original integrand is not.