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)$.