Count the trees
CJam (69 bytes)
]qi:X,{1e|,:):N{N\f{1$%!*}W$.*:+}%1$W%.*:+N,/+}/W\+_1,*X=\_W%.*:+-Y/z
Online demo
Explanation
The basic idea is to implement the generating function described in OEIS. The input \$0\$ is a nasty special case, but the final tweaks I made ended up producing \$-1\$ for that case, so the \$z\$ (for absolute value) tidies it up. That's the weirdest trick in here.
.*:+ is repeated three times, and looks like it could save a byte if extracted as {.*:+}:F~. However, this breaks with the special case \$0\$, because it doesn't execute the outer loop at all.
We use the auxiliary generating function for A000081, whose terms have the recurrence
a[0] = 0
a[1] = 1
For n >= 1, a[n+1] = (sum_{k=1}^n a[n-k+1] * sum_{d|k} d * a[d]) / n
I'm sure some languages have built-ins for the inverse Möbius transform \$\sum_{d\mid k} d \times a[d]\$, but CJam doesn't; the best approach I've found is to build an array mapping \$d\$ to k % d == 0 ? d : 0 and then do a pointwise multiplication with \$a\$ using .*. Note that here it is convenient to have built up \$a\$ starting at index 1, because we want to avoid division by zero when setting up the weights. Note also that if the two arrays supplied to the pointwise operation are not the same length then the values from the longer one are left untouched: therefore we have to either take the first \$k\$ terms of \$a\$ or make the array of weights go up to \$n\$. The latter seems shorter. So this inverse Möbius transform accounts for N\f{1$%!*}W$.*:+
If we call the result of the inverse Möbius transform M, we now have $$a[n+1] = \frac{1}{n} \sum_{k=1}^n a[n-k+1] \times M[k]$$
The numerator is obviously a term from a convolution, so we can handle it by reversing either a copy of \$a\$ or \$M\$ and then taking a pointwise multiplication and summing. Again, our index varies from \$1\$ to \$n\$, and in addition we want to pair up indices which sum to \$n + 1\$, so it's again convenient to index \$a\$ from 1 rather than 0. We've now accounted for
qi:X,{ ,:):N{N\f{1$%!*}W$.*:+}%1$W%.*:+N,/+}/
The point of the auxiliary generating function is given by the formula section of A000055:
G.f.: A(x) = 1 + T(x) - T^2(x)/2 + T(x^2)/2, where T(x) = x + x^2 + 2*x^3 + ... is the g.f. for A000081.
In terms of \$a\$, this means that the output we seek is $$[x = 0] + a[x] + \frac{1}{2}\left(a[x/2] - \sum_{i=0}^n a[i] \times a[n-i]\right)$$
where for \$a[x/2]\$ with odd \$x\$ we get 0. The shortest way I've found to do that is to inflate with zeroes (1,*) and then take the Xth element (X=).
Note that the sum is yet another convolution, but this time it's convenient to index from 0. The obvious solution is 0\+, but there is a nice little optimisation here. Since \$a[0] = 0\$, two terms of the convolution are guaranteed to be zero. We could take this as an opportunity to index from 1, but the special case \$X = 0\$ gets ugly. If we instead do W\+, the convolution gives us \$-2a[x] + \sum_{i=0}^n a[i] \times a[n-i]\$ and after subtraction and division by \$2\$ we've handled the \$a[x]\$ term outside the sum as well.
So we've explained
qi:X,{ ,:):N{N\f{1$%!*}W$.*:+}%1$W%.*:+N,/+}/W\+_1,*X=\_W%.*:+-Y/
The remaining details relate to the special case. I originally followed the recurrence more strictly by starting out with 1] and iterating from \$N = 1\$ with
1]qi:X,1>{ ... }/
The result when \$X = 0\$ is that we compute \$a\$ as [-1 1]: one term more than we need. The inflation and convolution wind up giving a result of \$0\$. (Perhaps better than we deserve - it's what we should have, because we haven't done anything about the term \$[x = 0]\$). So we fix it up for three chars either as a final X!+ or using the standard "fallback" technique as 1e|.
To get the "right" length of \$a\$ we need to avoid supplying that initial \$1\$, and instead to produce it from the main loop with \$N = 0\$. A completely straightforward
]qi:X,{ ... /+}/
obviously gives division by zero. But if we try
]qi:X,{1e| ... /+}/
then it works. We get
e# Stack: [] 0
1e| e# Stack: [] 1
,:):N e# Stack: [] [1]
{ e# We only execute this loop once
N\f{1$%!*} e# 1 divides 1, so stack: [] [1]
W$.* e# Remember: if the two arrays supplied to the pointwise operation
e# are not the same length then the values from the longer one are
e# left untouched. Stack: [] [1]
:+ e# Fold over a singleton. Stack: [] 1
}% e# And that was a map, so stack: [] [1]
1$W%.*:+ e# Another [1] [] .*:+, giving the same result: 1
N,/ e# 1 / 1 = 1
+ e# And we append 1 to a giving [1]
which produces exactly the value we require.
Now if \$X = 0\$ that main loop never executes, so after we hack in a \$-1\$ at the bottom we get [-1]. We end up computing \$(-1 - \frac12(-1 \times -1)) = -1\$, which isn't correct. But whereas correcting \$0\$ to \$1\$ took three chars, correcting \$-1\$ to \$1\$ only takes one: z gives the absolute value.