Counting with trees

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.

Another way to find the value of the sum $$ S:=\sum_{T} \prod_{i=1}^{n+1} d_i(T)! $$ using Cayley formula $$ \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} $$ is using not differentiating, as in the proof by Sam Hopkins and Ilya Bogdanov, but integrating: if you want to get $d!$ from $x^d$, you multiply by $e^{-x}$ and integrate from 0 to $\infty$. That said, we get $$S=\int_{[0,\infty)^{n+1}} x_1x_2\cdots x_{n+1} (x_1+x_2+\cdots+x_{n+1})^{n-1}e^{-x_1-x_2-\ldots-x_{n+1}}dx_1dx_2\ldots dx_{n+1}.$$ To find this integral, denote $$I_k=\int_{[0,\infty)^{n+1}} x_1x_2\cdots x_{n+1} (x_1+x_2+\cdots+x_{n+1})^{k}e^{-x_1-x_2-\ldots-x_{n+1}}dx_1dx_2\ldots dx_{n+1},$$ so, $S=I_{n-1}$. Using a natural change of variables $x_1+\ldots+x_{n+1}=y$, $x_1=t_1y$, $x_2=t_2y, \dots, x_n=t_ny$ we get $$I_k=\int_{0}^\infty dy\cdot e^{-y} y^{k+2n+1}\int_{t_i>0,t_1+\ldots+t_n<1}t_1\ldots t_n(1-t_1-\ldots-t_n)dt_1\ldots dt_n\\=(k+2n+1)!C,$$ where $C$ is some integral not depending on $k$. Since $I_0=(\int_0^\infty xe^{-x}dx)^{n+1}=1$, we get $C=1/(2n+1)!$ and therefore $$S=I_{n-1}=\frac{(3n)!}{(2n+1)!}=n! \, \frac{1}{2n+1}\binom{3n}{n}.$$

Well, and for what it worth an elementary proof. Let $\mathcal{F}_k$ denote the set of forests with $k$ trees, vertex set $\{1,2,\ldots,n+1\}$, and one labelled root in every component. Define the weight $w(F)$ of a forest $F\in \mathcal{F}_k$ as $\prod_{i=1}^{n+1} x_i!$, where $x_i=\deg_F(i)+\mathbf{1}(i\,\text{is a root in its component})$. Denote $s_k$ the sum of weights over $\mathcal{F}_k$. Note that $s_1=(3n+1)\times$ (your sum), since choosing a root in a fixed tree always gives a factor $\sum_i (d_i+1)=2n+(n+1)=3n+1$. I claim that $s_k=\frac{(3n-k+2)!}{(2n+1)!}{n\choose k-1}$ for all $k=1,2,\ldots,n+1$. This is obvious for $k=n+1$, so by (backwards) induction it suffices to verify $$s_k(n+1-k)=k(3n-k+2)s_{k+1}\tag{$\heartsuit$}$$ for all $k=1,2,\ldots,n$. We prove $(\heartsuit)$ by a double counting (of course) over the following graph, which is typical for such problems. For $F\in \mathcal{F}_k$, remove an edge, you get $k+1$ components, one of them, say $C$, without a root. Choose a root in $C$ as the only vertex of the removed edge which belongs to $C$ (another vertex of the removed edge belongs to a piece which already has a root). You get an element $\tilde{F}\in\mathcal{F}_{k+1}$, call $F$ and $\tilde{F}$ friends. Thus, each $F\in \mathcal{F}_{k}$ has $n+1-k$ friends (one friend for every edge). Thus, $$ s_k(n+1-k)=\sum_{F\in\mathcal{F}_{k},\tilde{F}\in \mathcal{F}_{k+1}\,\text{are friends}} w(F). $$ To prove $(\heartsuit)$ it remains to prove that for fixed $\tilde{F}\in \mathcal{F}_{k+1}$ the sum of $w(F)$ over all its friends $F\in \mathcal{F}_{k}$ equals $k(3n-k+2)w(\tilde{F})$. This is easy: if $C_1,\ldots,C_{k+1}$ are components of $\tilde{F}$, $r_1,\ldots,r_{k+1}$ are their roots, then to get a friend of $\tilde{F}$ you choose an ordered pair of indices $i\ne j$ and join $r_j$ to arbitrary vertex $v\in C_j$. The weight of this forest equals $w(\tilde{F})\cdot (x_v+1)$, where $x_v$ is defined as above (for the forest $\tilde{F}$). The sum of $x_v+1$ over $C_i$ is $2(|C_i|-1)+|C_i|+1=3|C_i|-1$, and this must be multiplied by $k$, since for fixed $i$ there exist $k$ ways to choose $j$ not equal to $i$. Totally we get a factor $k(3\sum |C_i|-(k+1))=k(3n-k+2)$.