Proof for symmetric polynomials with two variables X, Y:
Let $f(X,Y)$ be a symmetric polynomial with degree n. By definition, $f(X,Y)=f(Y,X)$.
Now, $$f(X,Y)= \sum_{0 \le i\le n,0\le j \le n,i+j\le n} a_{ij}X^iY^j$$
for some $a_{ij}$, and $a_{ij}=a_{ji} \;\forall{0\le i \le n,0\le j\le n}.\;\;\;\;(*)$
If $\alpha,\beta \in{\mathbb R}, c$ are some arbitrary constants, and $f(X,Y), g(X,Y)$ are two symmetric polynomials, then $\alpha f(X,Y)+\beta g(X,Y)+c$ is also a symmetric polynomial. This obviously extends to cases with multiple symmetric functions, i.e.
$$\sum_{j=1}^{p} \alpha_j f_j(X,Y)+c \; \text{is symmetric} \;\forall{\alpha_j \in{\mathbb R}, f_j \;\text{are all symmetric}} \;\;\;\;(**)$$
If n is even, then let $n=2m$
$$(X+Y)^n=X^n+\binom n 1 X^{n-1}Y+ \;...\;+\binom n mX^mY^m+ \;...\;+\binom n {n-1} XY^{n-1}+Y^n$$
If n is odd, then let $n=2m+1$
$$(X+Y)^n=X^n+\binom n 1 X^{n-1}Y+ \;...\;+\binom n mX^{m+1}Y^m+\binom n {m+1}X^mY^{m+1}+ \;...\;+\binom n {n-1} XY^{n-1}+Y^n$$
It is clear that $(X+Y)^n$ is symmetric $\forall{n\in {\mathbb N}}$. Moreover, $\forall 0\le k\le n,$ the $k^{th}$ term is the same as the $(n+2-k)^{th}$ term.
By $(*)$ and $(**)$ , it suffices to prove that every symmetric polynomial $f_0(X,Y)$ with at most 2 non-constant terms can be generated by $X+Y$ and $XY$, and that all polynomials with at most 2 non-constant terms generated by $X+Y$ and $XY$ are symmetric.
Suppose $\mu X^iY^j$ is symmetric, where $i\ge 0, j\ge 0$ are not both $0$. Then $X^iY^j=X^jY^i,$ which implies $i=j$.
Now, let us consider the polynomial $\;\eta_1 X^iY^j + \eta_2 X^rY^s$, where $i+j \ne r+s, i+j \ne 0, r+s \ne 0$. If it is symmetric, then $$\;\eta_1 X^iY^j + \eta_2 X^rY^s=\;\eta_1 X^jY^i + \eta_2 X^sY^r$$
This means either $X^iY^j=X^jY^i,X^rY^s=X^sY^r$ or $\eta_1=\eta_2, X^iY^j=X^sY^r, X^jY^i=X^rY^s$. Therefore, in the first case, $i=j, r=s$, and in the second case, $i=s,j=r$. Hence, all symmetric polynomials with two non-constant terms are of the form $\eta X^iY^i+\omega X^jY^j+c \;\text{or}\; \gamma (X^iY^j+X^jY^i)+c$.
Now, $\mu X^iY^i$ and $\eta X^iY^i+\omega X^jY^j$ can clearly be generated by $XY$, so we only have to check $\gamma (X^iY^j+X^jY^i)\;,$ where $i\ne j$.
Without loss of generality, suppose that $i\gt j$, so we have $$\gamma (X^iY^j+X^jY^i)=\gamma X^jY^j(X^{i-j}+Y^{i-j})$$. Now we only have to check if $$X^{i-j}+Y^{i-j}$$ can be generated by $X+Y$ and $XY$.
Let $i-j=q$. When $q=1$, $X^q+Y^q$ can trivially be generated by $X+Y$, since it's itself $X+Y$. When $q=2, X^q+Y^q=X^2+Y^2=(X+Y)^2-2XY,$ so $X^q+Y^q$ can be generated by $X+Y$ and $XY$ as well.
Notice that
$$X^{q+1}+Y^{q+1}=(X^q+Y^q)(X+Y)-(XY^q+YX^q)=(X^q+Y^q)(X+Y)-XY(X^{q-1}+Y^{q-1})$$
Therefore, $X^{q+1}+Y^{q+1}$ can be generated by $X+Y,XY,X^{q-1}+Y^{q-1},$ and $X^q+Y^q$. Hence, by induction,$\forall{q\in {\mathbb N}}, X^q+Y^q$ can be generated by $X+Y$ and $XY$. This means
$$\gamma(X^iY^j+X^jY^i)$$ can be generated by $X+Y$ and $XY$, hence all symmetric polynomials $f(X,Y)$ can be generated by $X+Y$ and $XY$.
All polynomials of degree n generated by $X+Y$ and $XY$ are of the form
$$\sum_{0 \le i'\le n,0\le j' \le n,i'+j'\le n} \kappa_{i'j'}(X+Y)^{i'}(XY)^{j'}$$
Each term of this sum is symmetric, so the sum is symmetric as well. Hence, every polynomial generated by $X+Y$ and $XY$ is symmetric, and every symmetric polynomial can be generated by $X+Y$ and $XY$.
If we make the further assumption that $f$ is homogeneous, then the solutions are the alternating polynomials, as I show below (the general case seems much harder). Note that those polynomials always have a symmetric square.
Suppose $f$ is an homogeneous solution, with total degree $d \gt 0$. Then, for any $k\geq 0$, $f^k$ is also homogeneous with total degree $kd$. If we write
$Q=\sum_{j=0}^d q_jx^j$ and $q_d\neq 0$, we see that the components $q_jf^j$ do not overlap, so that $f$ is symmetric iff each $q_jf^j$ is.
In particular, $q_df^d$ is symmetric and we may assume $Q=x^d$.
So we have the identity $f(x_1,x_2,\ldots,x_n)^d = f(x_{\sigma(1)},x_{\sigma(2)},\ldots,x_{\sigma(n)})^d$ in ${\mathbb Q}[x_1,x_2,\ldots,x_n]$ for every $\sigma\in {\mathfrak S}_n$, whence
$f(r_{\sigma(1)},\ldots,r_{\sigma(n)})=\pm f(r_1,r_2,\ldots,r_n)$ for any $(r_1,r_2,\ldots,r_n)\in {\mathbb Q}^n$. So there is a map $\lambda:{\mathfrak S}_n \times {\mathbb Q}^n \to \lbrace \pm 1 \rbrace$ such that
$$
f(r_{\sigma(1)},\ldots,r_{\sigma(n)})=\lambda(\sigma,r_1,r_2,\ldots,r_n) f(r_1,r_2,\ldots,r_n) \tag{1}
$$
Now, fix $\sigma\in {\mathfrak S}_n$ and $s=(r_1,\ldots,r_{n-1})\in {\mathbb Q}^{n-1}$, and consider the partial map $\lambda_{\sigma,s}:{\mathbb Q} \to \lbrace \pm 1 \rbrace$ defined by $\lambda_{\sigma,s}(r_n)=\lambda(r_1,\ldots,r_{n-1},r_n,\sigma)$. Since the image of $\lambda_{\sigma,s}$ is finite and $\mathbb Q$ is infinite, there is a value $\mu(\sigma,s)\in \lbrace \pm 1 \rbrace$ that is attained infinitely often. The identity $f(s.\sigma,r_{\sigma(n)})=\mu(\sigma,s)f(s,r_n)$ then holds for infinitely many $r_n$, and hence holds for all $r_n$ since $f$ is a polynomial.
Repeating this argument and using induction on $n$, we eventually see that we can take $\lambda$ to be dependent on $\sigma$ only and not on $r_1,\ldots,r_n$ :
$$
f(r_{\sigma(1)},\ldots,r_{\sigma(n)})=\lambda(\sigma) f(r_1,r_2,\ldots,r_n) \tag{2}
$$
Since we assume $f\neq 0$, there is a $(\rho_1,\rho_2,\ldots,\rho_n)$ with $f(\rho_1,\rho_2,\ldots,\rho_n)\neq 0$. Then we can write :
$$
\lambda(\sigma)=\frac{f(\rho_{\sigma(1)},\rho_{\sigma(2)},\ldots,\rho_{\sigma(n)})}{f(\rho_1,\rho_2,\ldots,\rho_n)} \tag{3}
$$
and it easily follows that $\lambda$ is a group homomorphism ${\mathfrak S}_n\to \lbrace \pm 1 \rbrace$. So, $\lambda$ either the identity or the signature homomorphism, and $f$ is an alternating polynomial.
Best Answer
I suppose I should actually come back to answer this. Main idea: for $1\le d<n$ we have
$$\Bbb Z[x_1,\dots,x_n]^{S_d}=\bigoplus_{j=0}^d \Bbb Z[x_1,\dots,x_n]^{S_{d+1}}x_{d+1}^j. \tag{$\circ$}$$
They're both equal to $\Bbb Z[x_1,\dots,x_n]^{S_{d+1}}[x_{d+1}]$. Using $(\circ)$ we can begin unpeeling $\Bbb Z[x_1,\dots,x_n]$ by setting $d=1,2,\dots,n-1$. What I found interesting about this argument is that the outer layers of the onion are the smaller indices rather than the larger - initially I believed it'd be the reverse.