The answer is YES to all your questions : yes it is already known, and I'm pretty certain that I've already seen here on MSE if not a proof, at least a mention of it. But I can't seem to find it back so I rewrite a proof here.
Let $(f_k)_{k\geq 0}$ be a sequence of polynomials with $\deg(f_k)=k$
and $f_i(f_j(x))=f_j(f_i(x))$ for any $i,j$. Using a suitable conjugation, we may assume
$$f_0(x)=1\tag{1}$$
Then, for any $i$ we have $f_i(1)=1$, or in other words
$$f_i(x)=1+(x-1)g_i(x)\tag{2}$$
for some $g_i$ of degree $i-1$. There are two constants $c$ and $d$ such that $g_2(x)=cx+d$ (note that $c\neq 0$). It follows that $f_2(x)=cx^2+(d-c)x+1-d$, and hence $f_1(f_2(x))=g_1f_2(x)+1-g_1=(cg_1)x^2+g_1(d-c)x+1-g_1d$ while
$f_2(f_1(x))=(cg_1^2)x^2+g_1(d+c-2g_1)x+(g_1^2-g_1)c+1-g_1d$. Comparing the leading coefficients and using $cg_1\neq 0$, we deduce that $g_1=1$, so
$$
f_1(x)=x \tag{3}
$$
There are three constants $k,l,m$ such that $g_3=kx^2+lx+m$. It follows that $f_3(x)=kx^3 + (l-k)x^2 + (m-l)x + 1-m$, and hence
$$
\begin{array}{lcrrrr}
kf_2(x)^3 & = & (kc^3)x^6 + 3kc^2(d-c)x^5 & + 3kc(c^2+d^2-3cd+c)x^4
& + k(d-c)(c^2+d^2-8cd+6c)x^3 & + \ldots \\
(l-k)f_2(x)^2 & = & & + (l-k)c^2 x^4
& + 2c(l-k)(d-c)x^3 & + \ldots \\
f_3(f_2(x)) & = & (kc^3)x^6 + 3kc^2(d-c)x^5 & + c\big(lc+k(3c^2+3d^2-9cd+2c)\big) x^4
& + (d-c)\big(2cl+k(c^2+d^2-8cd+4c)\big)x^3 & + \ldots \\
\end{array}
$$
Similarly,
$$
\begin{array}{lcrrr}
cf_3(x)^2 & = & (k^2c)x^6 + 2kc(l-k)x^5 + c(k^2+l^2+2mk-4kl)x^4
& + 2c(-l^2+(k+m)l+k(1-2m))x^3 & + \ldots \\
(d-c)f_3(x) & = &
& + k(d-c)x^3 & + \ldots \\
f_2(f_3(x)) & = & (k^2c)x^6 + 2kc(l-k)x^5 + c(k^2+l^2+2mk-4kl)x^4
& + \big(kd+2c(-l^2+(k+m)l+k(\frac{1}{2}-2m))\big)x^3 & + \ldots \\
\end{array}
$$
We thus obtain the system
$$
\left\lbrace
\begin{array}{lcl}
kc^3 &=& k^2c \\
3kc^2(d-c) &=& 2kc(l-k) \\
c\big(lc+k(3c^2+3d^2-9cd+2c)\big) &=& c(k^2+l^2+2mk-4kl) \\
(d-c)\big(2cl+k(c^2+d^2-8cd+4c)\big) &=& kd+2c(-l^2+(k+m)l+k(\frac{1}{2}-2m)) \\
\end{array}\right.
$$
Remembering that $ck\neq 0$, the first line yields $k=c^2$, the second then yields $l=\frac{c}{2}(3d-c)$, and the third yields $m=\frac{-c^2+6(1-d)c+3d(d+2)}{8}$. The fourth equation then becomes equivalent to $c^5 + (3d - 6)c^4 + (3d^2 - 12d + 8)c^3 + (d^3 - 6d^2 + 8d)c^2=0$, which nicely factorizes as $c^2(d+c-4)(d+c-2)(d+c)=0$. So $d$ must be one of $4-c,2-c$ or $-c$.
When $d=4-c$, putting $\sigma(x)=\frac{(4-d)x-(2-d)}{2}$, we have
$\sigma f_2\sigma^{-1}(x)=2x^2-1$ and $\sigma f_3\sigma^{-1}(x)=4x^3-3x$.
When $d=2-c$, putting $\sigma(x)=(2-d)x-(1-d)$, we have
$\sigma f_2\sigma^{-1}(x)=x^2$ and $\sigma f_3\sigma^{-1}(x)=x^3$.
When $d=-c$, putting $\sigma(x)=d(1-x)$, we have
$\sigma f_2\sigma^{-1}(x)=x^2$ and $\sigma f_3\sigma^{-1}(x)=x^3$.
So we have shown that (up to conjugation) we can assume either (a) $f_2(x)=x^2, f_3(x)=x^3$, (b) $f_2(x)=T_2(x)$, $f_3(x)=T_3(x)$.
Suppose we are in case (a). Let $n\geq 4$. Then $f_n(x^2)=f_n(x)^2$. Comparing the leading coefficients, we see that $f_n$ is monic. Let $z\in{\mathbb C}$ be a complex root of $f_n$. Since $f_n(x^2)=f_n(x)^2$, we see that $z^2$ is a root of $f_n$ also, and by induction every $z^{2^m}$ for $m\geq 1$ is a root of $f_n$. Since there are only finitely many roots, by the pigeonhole principle we can find $m\lt m'$ such that $z^{2^m}=z^{2^{m'}}$, and if $z\neq 0$ then $z$ must be a root of unity. But a root of unity is of the form $exp(2pi r)$ where $r \in {\mathbb Q}/{\mathbb Z}$. Denote by $R$ the (possibly empty) subset of ${\mathbb Q}/{\mathbb Z}$ of all arguments that can be obtained in this way. Because of $f_n(x^2)=f_n(x)^2$, $R$ is stable by division by $2$, so it is either empty or infinite (which is impossible). Thus $0$ is the only complex root of $f_n$, $f_n(x)=x^n$ and we are done in this case.
Suppose we are in case (b). Let $n\geq 4$. Then $f_n(T_2(x))=T_2(f_n(x))$. Comparing the leading coefficients, we see that the leading coefficient of $f_n$ is $2^{n-1}$. Write $f_n(x)=\sum_{k=0}^n a_kx^{n-k}$ (so $a_0=2^{n-1}$). Then, for $1\leq d \leq n$ in the expansion of $f_n(T_2(x))-T_2(f_n(x))$, the coefficient of degree $2n-d$ in $x$ is of the form $-4a_0a_d+h(a_0,a_1,\ldots,a_{d-1})$ where $h$ is a polynomial in $a_0,a_1,\ldots,a_{d-1}$. Thus the value of $a_d$ can be deduced from the values of $a_0,a_1,\ldots,a_{d-1}$. It follows that there is at most one polynomial of degree $n$ commuting with $T_2$. But we already know that $T_n$ has this property, so $f_n=T_n$ and we are done.