A conjecture relating Multiple Zeta Values and the Polya Enumeration Theorem Let me state my motivation. I believe that the
Polya Enumeration Theorem
and Multiple Zeta Values (the classic  being the  Basel problem and  the values of  the Riemann
zeta  function at  the even  integers) are  among the  most intriguing
objects in mathematics.

With this message I present  a conjecture that connects these two (PET
and MZVs).  Most likely there  is a proof somewhere in the literature,
which I invite  readers to submit in text form if  it admits a compact
presentation or as a reference if it does not.

To understand this  conjecture you need to learn  about multiset cycle
indices, which I presented in a different context at this
MSE link.
Multisets and their  cycle indices form a combinatorial  species like  any other,
e.g.  cycles, sequences  and sets.  They represent  multisets  and are
identified by the partition  that corresponds to the multiplicities of
the  elements  of  the  multiset,  e.g.  $\mathfrak{M}_{1,2,3}$  is  a
multiset that contains three elements,  one of which in two copies and
another in three.  Substitution of an OGF into  a multiset cycle index
yields the generating function of the multiset.
There   is  a   special  multiset   cycle  index   which   is  written
$\mathfrak{M}_{1,1,1,\ldots,1}$  which gives  the  multiset where  all
elements are  unique, i.e. the species  of sets. This  cycle index has
been known for quite some time and is used to compute the set operator
$\mathfrak{P}.$  The   cycle  index  $Z(P_n)$  of   the  set  operator
$\mathfrak{P}_{=n}$ is  the difference between the cycle  index $Z(A_n)$ of
the alternating  group and the  cycle index $Z(S_n)$ of  the symmetric
group on $n$ elements. It admits the following simple recursive definition:
$$Z(P_0) = 1 \quad\text{and}\quad
Z(P_n) = \frac{1}{n}
\sum_{l=1}^n (-1)^{l+1} a_l Z(P_{n-l}).$$
Here are the cycle indices $Z(P_3), Z(P_4)$ and $Z(P_5):$
$$\begin{array}{|l|l|}
\hline
Z(P_3) &
\frac{1}{6}\,{a_{{1}}}^{3}-1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}}\\
\hline
Z(P_4) &
\frac{1}{24}\,{a_{{1}}}^{4}-1/4\,a_{{2}}{a_{{1}}}^{2}
+1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}-1/4\,a_{{
4}}\\
\hline
Z(P_5) &
{\frac {1}{120}}\,{a_{{1}}}^{5}-\frac{1}{12}\,a_{{2}}{a_{{1}}}^{3}
+1/6\,a_{{3}}{a_{{1}}}^{2}+1/8\,a_{{
1}}{a_{{2}}}^{2}-1/4\,a_{{4}}a_{{1}}
-1/6\,a_{{2}}a_{{3}}+1/5\,a_{{5}}\\
\hline\end{array}$$
With these definitions we are  now ready to state the conjecture which
is quite simply that
$$\large \color{blue}{\zeta(s, s, s, \ldots, s) = 
Z(P_n)(\zeta(s), \zeta(2s), \ldots, \zeta(ns))}$$
i.e.   the MZV  of a  unique argument  $s$ repeated  $n$ times  is equal to the
substituted  cycle  index  $Z(P_n)$  of the set operator $\mathfrak{P}_{=n}$ according  to  the  rule  $a_l  =
\zeta(ls).$
Note that the case $n=3$ is  found on the Wikipedia page for MZVs that
I linked to in the introduction.

Example. The conjecture says e.g. that
$$\zeta(5,5,5,5) =
1/24\,  \zeta \left( 5 \right) ^{4}-1/4\,\zeta \left( 10 \right)
 \zeta \left( 5 \right)  ^{2}+1/3\,\zeta \left( 15 \right) \zeta
 \left( 5 \right) +1/8\,  \zeta \left( 10 \right) ^{2}-1/4\,\zeta
 \left( 20 \right).$$

Important observation.  The reader may well  ponder the statement
of the conjecture and say  that there is nothing to  prove here, which I  will accept as
proof, if accompanied by a brief explanation why this should be so.
Addendum. An interesting application of the above result is found at this MSE link.
 A: This actually has nothing to do with zeta functions, it is pure symmetry - a special case of
$$e_n=Z(P_n)(p_1,p_2,\cdots,p_n) \tag{1}$$
where $e_k$ and $p_k$ are the elementary and power-sum symmetric polynomials, in infinitely many variables $x_1,x_2,\cdots$, respectively. For zeta values we simply evaluate at $x_k:=k^{-s}$.
The recursive identity you wrote is the content of Newton's identities, which yields explicitly
$$e_n=\frac{1}{n}\sum_{l=1}^n(-1)^{l+1}p_le_{n-l}. \tag{2}$$
A more explicit formula is possible using generating functions. Observe
$$\begin{align} \sum_{n=0}^\infty T^ne_n & = \prod_{i=1}^\infty(1+Tx_i) \tag{3} \\
& = \exp\left[\sum_{i=1}^\infty \log(1+Tx_i)\right] \tag{4} \\
& = \exp\left[\sum_{i=1}^\infty \sum_{\ell=1}^\infty (-1)^{\ell+1}\frac{(Tx_i)^\ell}{\ell}\right] \tag{5} \\
& = \exp\left[\sum_{\ell=1}^\infty(-1)^{\ell+1}\frac{T^\ell}{\ell}\sum_{i=1}^\infty x_i^\ell\right] \tag{6} \\
& = \exp\left[\sum_{\ell=1}^\infty(-1)^{\ell+1}\frac{T^\ell}{\ell}p_\ell\right] \tag{7} \\
& = \prod_{\ell=1}^\infty\exp\left[(-1)^{\ell+1}\frac{T^\ell}{\ell}p_\ell\right] \tag{8} \\
& = \prod_{\ell=1}^\infty \left[\sum_{c_\ell=0}^\infty \frac{1}{c_\ell!}\left((-1)^{\ell+1}\frac{T^\ell}{\ell}p_\ell\right)^{c_\ell}\right] \tag{9} \\
& = \sum_{n=0}^\infty T^n \sum_{1c_1+2c_2+\cdots=n} \prod_{\ell=1}^\infty \frac{(-1)^{(\ell+1)c_\ell}}{c_\ell!}\left(\frac{p_\ell}{\ell}\right)^{c_\ell} \tag{10} \\
& = \sum_{n=0}^\infty \frac{T^n}{n!} \sum_{1c_1+\cdots+nc_n=n} (-1)^{c_2+c_4+\cdots}\frac{n!}{1^{c_1}2^{c_2}\cdots n^{c_n}c_1!c_2!\cdots c_n!}p_1^{c_1}\cdots p_n^{c_n} ~~ \tag{11} \\ 
& = \sum_{n=0}^\infty \frac{T^n}{n!}\sum_{\sigma\in S_n} {\rm sgn}(\sigma)p_1^{c_1(\sigma)}\cdots p_n^{c_n(\sigma)} \tag{12} \end{align}$$
where $c_k(\sigma)$ is the number of length $k$ cycles in $\sigma$'s disjoint cycle decomposition, the sign of $\sigma$ is equal to $(-1)^{c_2+c_4+\cdots}=(-1)^{2c_1+3c_2+4c_3+\cdots}$, and the number of permutations in $S_n$ with the cycle type $(\underbrace{n,\cdots,n}_{c_n},\cdots,\underbrace{1,\cdots,1}_{c_1})$ is given $n!/(1^{c_1}2^{c_2}\cdots n^{c_n}c_1!\cdots c_n!)$ (see here).
Therefore we may conclude
$$e_n=\frac{1}{n!}\sum_{\sigma\in S_n}{\rm sgn}(\sigma)p_1^{c_1(\sigma)}\cdots p_n^{c_n(\sigma)}. \tag{13}$$
The more general polynomials in your additional MZV identities are called the monomial symmetric functions $m_\mu$ for integer partitions $\mu$. These form a linear basis for the space of symmetric functions. I am not aware of a direct formula for monomial symmetric functions in terms of power sums, but in terms of standard constants and functions there is an indirect way: writing monomials in terms of Schur polynomials $s_\lambda$, and writing those with power-sums.
The Kostka numbers are defined by the relation $s_\lambda=\sum_\mu K_{\lambda\mu}m_\mu$, where the sum is over integer partitions. If one inverts the Kostka matrix we can write $m_\mu=\sum_\lambda J_{\mu\lambda}s_\lambda$ for some values $J_{\mu\lambda}$.
Schur polynomials in turn can be decomposed as
$$s_\lambda=\frac{1}{n!}\sum_{\sigma\in S_n}\chi_\lambda(\sigma)p_1^{c_1(\sigma)}\cdots p_n^{c_n(\sigma)} \tag{14}$$
where $\chi_\lambda$ is the irreducible character of $S_n$ corresponding to $\lambda$.
Another route for envelope calculation is to decompose $m_\mu$ into the variables $e_k$ (which is possible thanks to the fundamental theorem of symmetric polynomials) using this recursive algorithm, and then writing the $e_k$s in terms of $p_k$s as in $(13)$. This should be faster, although likely there are even faster, more direct ways of writing monomials in terms of power-sums. For more literature related to this you'll want to read deeper into the combinatorics of symmetric polynomials as well as the representation theory of symmetric groups.
Out of theoretical interest, representation theory illustrates $(14)$ is a trace formula:
$$\chi_{\Bbb S_\lambda(V)}(g)=\frac{1}{n!}\sum_{\sigma\in S_n} \chi_\lambda(\sigma)\chi_V(g^1)^{c_1(\sigma)}\cdots\chi_V(g^n)^{c_n(\sigma)}$$ where $\Bbb S_\lambda(V)=V^{\otimes n}\otimes_{\Bbb C[S_n]}M_\lambda$ is the Schur functor applied to $V$ considered as a representation of the group ${\rm GL}(V)$, and $M_\lambda$ is the irreducible representation of $S_n$ associated to $\lambda$.
