Denote $x = (x_1,...,x_n)$. I'm trying to prove the following:
$$\int_{S^{n-1}}x_1^2dS =\int_{S^{n-1}}x_k^2dS \; , \; 2\leq k\leq n $$
Intuitively this equality is due to the symmetry of the sphere, but I'm looking for a formal explanation. I thought about using the definition with a parametrization, but I'm not sure how to find a good paramterization of the sphere for that purpose.
A formal explanation uses the invariance of the measure $dS$ under orthogonal transformations. That is, given a continuous function $f:\mathbb{S}^{n-1}\to\mathbb{R}$ and an orthogonal matrix $U$ (that is, $U^tU=I$), one has $$\int_{\mathbb{S}^{n-1}}f(Ux)\,dS(x)=\int_{\mathbb{S}^{n-1}}f(x)\,dS(x)\quad (*)$$
In your case, take $f(x)=\langle x,e_1\rangle^2$, where $\langle\cdot\rangle$ denotes the standard inner product. Then $f(x)$ returns the square of the first coordinate of $x$. Pick any $1\leq k\leq n$. Take an orthogonal matrix $U$ such that $U^te_1=e_k$. (That is, rotate the vector $e_k$ to the vector $e_1$). Then, by invariance, you have $$ \begin{aligned} \int_{\mathbb{S}^{n-1}}x_1^2\,dS&=\int_{\mathbb{S}^{n-1}}\langle x,e_1\rangle^2\,dS\\ &=\int_{\mathbb{S}^{n-1}}\langle Ux,e_1\rangle^2\,dS =\int_{\mathbb{S}^{n-1}}\langle x,U^te_1\rangle^2\,dS\\ &=\int_{\mathbb{S}^{n-1}}x_k^2\,dS \end{aligned} $$