Given $Z_1, \ldots, Z_{n + 1} \text{ i.i.d.} \sim \text{Exp}(1)$, it can be shown that (most easily by characteristic function or moment generating function) $A_n \sim \text{Gamma}(i, 1), B_n \sim \text{Gamma}(n - i + 1, 1)$. Therefore, what to be proved is a corollary of the result below:
If $X \sim \text{Gamma}(k_1, \theta), Y \sim \text{Gamma}(k_2, \theta)$, $X$ and $Y$ are independent, then
$$\frac{X}{X + Y} \sim \text{Beta}(k_1, k_2). \tag{1}$$
The standard way of showing $(1)$ is evaluating the probability $F(t) = P[X/(X + Y) \leq t]$ for $0 < t \leq 1$ directly:
\begin{align}
F(t) &= P[X/(X + Y) \leq t] = \int_0^\infty P\left[\frac{x}{x + Y} \leq t\right]f_X(x)~\mathrm dx \\
&= \int_0^\infty P[Y \geq (t^{-1} - 1)x]f_X(x)~\mathrm dx \\
&= \int_0^\infty\left[\int_{(t^{-1} - 1)x}^\infty f_Y(y)dy\right] f_X(x)~\mathrm dx,
\end{align}
where $f_X, f_Y$ are densities of $X$ and $Y$ respectively. The density $f(t)$, of $X/(X + Y)$, can thus be obtained by taking derivative of $F$ with respect with $t$ under the integral as follows:
\begin{align}
& f(t) = F'(t) = \int_0^\infty f_Y((t^{-1} - 1)x)\frac{x}{t^2}f_X(x) ~\mathrm dx \\
=& \int_0^\infty \frac{x}{t^2}\frac{1}{\Gamma(k_2)\theta^{k_2}}(t^{-1} - 1)^{k_2 - 1}x^{k_2 - 1}e^{-\frac{(t^{-1} - 1)x}{\theta}} \times
\frac{1}{\Gamma(k_1)\theta^{k_1}}x^{k_1 - 1}e^{-\frac{x}{\theta}}~\mathrm dx \\
=&\frac{\Gamma(k_1 + k_2)(1 - t)^{k_2 - 1}t^{k_1 + k_2}}{\Gamma(k_1)\Gamma(k_2)t^{k_2 + 1}}
\boxed{\int_0^\infty\frac{1}{\Gamma(k_1 + k_2)(\theta t)^{k_1 + k_2}}x^{k_1 + k_2 - 1}e^{-\frac{x}{\theta t}}~\mathrm dx} \\
=& \frac{(1 - t)^{k_2 - 1}t^{k_1 - 1}}{B(k_1, k_2)},
\end{align}
which coincides the density function of a $\text{Beta}(k_1, k_2)$ random variable. The integral inside the box is unity because it is the integral of the density function of a $\text{Gamma}(k_1 + k_2, \theta t)$ random variable. This completes the proof.
