The notation is potentially confusing, but the end result is correct. Essentially this claim is the relationship between the area and the energy of of a conformal map. Here is a slightly expanded proof of their computation:
Let $X_i:\mathbb{R}^{n+1}\to\mathbb{R}$ denote the $i$-th coordinate function.
Suppose that $\phi:M^2\to \mathbb{S}^n$ is a conformal map from a surface. Then, they argue that there is a new conformal map $F:=g\circ \phi:M\to \mathbb{S}^n$ which balances the $X_i$ in the sense that $$ \int_M X_i\circ F = 0 $$ for each $i$. Thus, $X_i\circ F$ is an acceptable test function for the variational characterization of the first eigenfunction.
Choose isothermal coordinates, $y_i$ near some point $p\in M$. That is, $g = \lambda^2 \delta_{\mathbb{R}^2}$. Because $F$ is conformal and because we have chosen isothermal coordinates, we see that the vectors in $\mathbb{R}^{n+1}$ $$ F_*\left(\frac{\partial}{\partial y_i}\right) $$ for $i=1,2$ are orthogonal and of the same length $\lambda^{-1}$ (i.e. not necessarily unit vectors). Then, we may compute the area of the image of $F$ as (lets assume that the coordinates $y_i$ cover $M$, otherwise this computation is valid locally and then must be patched together) \begin{align*} \mathrm{Area}(F)&=\int_M \Vert F_*(\partial/\partial y_1)\wedge F_*(\partial/\partial y_2)\Vert dy_1\wedge dy_2\\ & = \int_M \sqrt{ \Vert F_*(\partial/\partial y_1) \Vert^2 \Vert F_*(\partial/\partial y_2) \Vert^2 - \langle F_*(\partial/\partial y_1), F_*(\partial/\partial y_2)\rangle} dy_1 \wedge dy_2\\ & = \int_M { \Vert F_*(\partial/\partial y_1) \Vert \Vert F_*(\partial/\partial y_2) \Vert } dy_1 \wedge dy_2\\ & = \frac 12 \sum_{i=1}^2\int_M \Vert F_*(\partial/\partial y_i) \Vert^2 dy_1 \wedge dy_2\\ & = \frac 12 \sum_{i=1}^2 \sum_{j=1}^{n+1} \int_M\left( \frac{\partial F_j}{\partial y_i}\right)^2 dy_1 \wedge dy_2\\ & = \frac 12 \sum_{j=1}^{n+1} \int_M \lambda^{-2}\Vert \nabla (X_j\circ F)\Vert^2 dy_1 \wedge dy_2\\ & = \frac 12 \sum_{j=1}^{n+1} \int_M \Vert \nabla (X_j\circ F)\Vert^2 dV_g\\ \end{align*}
Let me explain these steps: The first equality is by definition, as is the second. Then, the third follows because the two vectors are orthogonal, so their inner product vanishes. Furthermore, they have the same length, so the fourth equality holds. The next is basically by definition of the Euclidean norm. Then, we change this into the $g$-gradient of $X_j\circ F$ by definition of the gradient and by the choice of the isothermal coordinates. Finally, we use the expression for the volume form in coordinates.
