Now that we have the machinery of species, the plane trees problem cracks right open. We will use the following species in our computation:

• The species $R$ of rooted plane trees, that is, plane trees with one vertex singled out as the root. Notice that there are $n$ times as many rooted plane trees as plane trees on $n$ vertices, so $$\widetilde{R}(x)=\sum_{n\ge 0} \frac{np_n}{n!}x^n=\sum_{n\ge 0} \frac{p_n}{(n-1)!}x^n.$$
• The species $T$ of ordered rooted trees, that is, a rooted tree such that the children of each node starting from the root have a specified linear ordering.
• The species $X$ defined by $X([n])=\begin{cases}1 & n=1 \\ \emptyset & n\neq 1\end{cases}$. Notice that $$\widetilde{X}(x)=x.$$
• The species $C$ of cyclic orderings of $[n]$. Notice that $$\widetilde{C}(x)=\sum_{n\ge 0}\frac{1}{n}x^n=-\ln(1-x).$$
• The species $L$ of linear orderings of $[n]$. Notice that $$\widetilde{L}(x)=\frac{1}{1-x}.$$

Now, notice that by partitioning a rooted plane tree into its root and its collection of branches with a cyclic ordering on the branches, we have $R=X\cdot (C\odot T)$, as each of the branches can be thought of as a rooted tree along with a linear ordering the children of each node. An example with $1$ as the root, and its branches circled, is shown below.

Taking exponential generating functions, we have $$\begin{equation}\widetilde{R}(x)=x\cdot (-\ln(1-\widetilde{T}(x))),\label{eq:one}\end{equation}$$ so we wish to find $\widetilde{T}(x)$. By a similar species argument, $T=X\cdot (L\odot T)$, and so $$\widetilde{T}(x)=\frac{x}{1-\widetilde{T}(x)}.$$ Solving, we find $$\widetilde{T}(x)=\frac{1-\sqrt{1-4x}}{2},$$ and combining this with equation $\eqref{eq:one}$ we find $$\widetilde{R}(x)=-x\ln\left(\frac{1+\sqrt{1-4x}}{2}\right).$$

We’re almost there: now we just have to extract the coefficients. We have $$\left(\frac{\widetilde{R}(x)}{x}\right)’=\sum_{n\ge 0} \frac{p_{n+2}}{n!} x^n,$$ and by our formula for $\widetilde{R}(x)$,
$$\begin{eqnarray*}
\left(\frac{\widetilde{R}(x)}{x}\right)’ &=& \frac{2}{1+\sqrt{1-4x}}\cdot\frac{1}{\sqrt{1-4x}} \\
&=& \frac{2}{\sqrt{1-4x}}\left(\frac{1-\sqrt{1-4x}}{4x}\right) \\
&=& \frac{1}{2x}(-1+(1-4x)^{-1/2}) \\
&=& \frac{1}{2x}\left(\sum_{k=1}^\infty (-1)^k\binom{-1/2}{k}4^k x^k\right) \\
&=& \sum_{k=0}^\infty (-1)^{k+1} 2^{2k+1} \binom{-1/2}{k} x^k \\
&=& \sum_{k=0}^\infty 2^{2k+1} \frac{(1/2)(3/2)\cdots ((2k+1)/2)}{(k+1)!} x^k \\
&=& \sum_{k=0}^\infty 2^k\frac{(2k+1)!!}{k+1} \frac{x^k}{k!} \\
\end{eqnarray*}$$ (Here $(2k+1)!!$ means $1\cdot 3 \cdot 5 \cdot\cdots\cdot (2k+1)$.)

And so, comparing coefficients, we conclude that $$p_{n+2}=2^n\frac{(2n+1)!!}{(n+1)}$$ for all $n\ge 0$. Indeed, this checks out with our previous computations of $p_2=1$, $p_3=3$, $p_4=20$, and $p_5=210$.

So species have saved the day… but now that we know the formula for $p_n$, can we find a direct combinatorial proof? I haven’t found one yet, so feel free to comment with ideas!