Here's a start. Your claim is equivalent to (and easier to understand as)
$$ \sum_{T} \prod_{i=1}^{n+1} d_i(T)! = n! \, \frac{1}{2n+1}\binom{3n}{n}$$
where the sum is over all labeled trees $T$ on vertex set $[n+1]:=\{1,2,\ldots,n+1\}$ and $d_i(T)$ is the degree of vertex $i$ in $T$.

Recall that the (generalized) Cayley formula says that
$$ \sum_{T} x_1^{d_1(T)}x_2^{d_2(T)}\cdots x_{n+1}^{d_{n+1}(T)} = x_1x_2\cdots x_{n+1} (x_1+x_2+\cdots+x_{n+1})^{n-1}$$
where again the sum is over all labeled trees $T$ on $[n+1]$ and $d_i(T)$ is as before.

Hence, defining the differential operator $\Psi$ by
$$ \Psi = \sum_{(\alpha_1,\ldots,\alpha_{n+1})\vDash 2n} \frac{\partial^{2n}}{\partial x_1^{\alpha_1} \cdots \partial x_{n+1}^{\alpha_{n+1}}}$$
where the sum is over all compositions $(\alpha_1,\ldots,\alpha_{n+1})$ of $2n$ satisfying $1 \leq \alpha_i \leq n$ for all $i$, your claim becomes
$$\Psi(x_1\cdots x_{n+1} (x_1+\cdots+x_{n+1})^{n-1}) = n!\,\frac{1}{2n+1}\binom{3n}{n}.$$
This last expression at least avoids any discussion of trees.