With the question by the OP asking for a reference I will try to do
just that, providing links to the OEIS. Surprisingly enough even the
OEIS does not offer the usual variety of references here, suggesting
that this problem is open. We compute generating functions $T_{\le
h}(z)$ for the height being at most $h$ and the desired count is then
given by $T_{\le h}(z)-T_{\le h-1}(z).$
We have for the labeled the combinatorial class
$$\mathcal{T}_{\le h} =
\mathcal{Z}\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}\times \textsc{SET}(\mathcal{T}_{\le h-1}).$$
For the labeled case this translates into
$$T_{\le h}(z) = z\exp T_{\le h-1}(z)$$
with $T_{\le 0}(z) = 0$ and $T_{\le 1}(z) = z.$
For example with $h=3$ we get
$$T_{\le 3}(z) = z \exp(z \exp(z))$$
which yields the sequence OEIS A052512:
$$1, 2, 9, 40, 205, 1176, 7399, 50576, \ldots$$
and for $h=4$ we get
$$T_{\le 4}(z) = z(\exp( z \exp(z \exp(z)))$$
which yields the sequence OEIS A052513:
$$1, 2, 9, 64, 505, 4536, 46249, 526352, \ldots$$
We have for the unlabeled the combinatorial class
$$\mathcal{T}_{\le h} =
\mathcal{Z}\times \textsc{MSET}(\mathcal{T}_{\le h-1}).$$
which translates into generating functions, as follows,
$$T_{\le h}(z) = z\exp
\left(\sum_{l\ge 1}
\frac{T_{\le h-1}(z^l)}{l}\right).$$
where $T_{\le 0}(z) = 0$ and $T_{\le 1}(z) = z$ as before.
We get for
$$T_{\le 2}(z) =
z\exp \left(\sum_{l\ge 1} \frac{z^l}{l}\right)
= z\exp\log\frac{1}{1-z} = \frac{z}{1-z}.$$
These are root nodes with singletons attached to them.
Continuing with the recursion we obtain
$$T_{\le 3}(z) = z\exp
\left(\sum_{l\ge 1}
\frac{1}{l}\sum_{q\ge 1} z^{ql}\right)
= z\exp
\left(\sum_{q\ge 1}
\sum_{l\ge 1} \frac{1}{l} z^{ql}\right)
\\ = z\exp\sum_{q\ge 1}\log\frac{1}{1-z^q}
= z \prod_{q\ge 1} \frac{1}{1-z^q}.$$
which yields the sequence OEIS A000041:
$$1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, \ldots$$
(partition numbers, not surprisingly, as we allocate one node for the
root and the rest is a partition into singleton fans of height at most
two, of which there is only one per node count).
From $T_{\le 4}(z)$ we obtain the sequence
OEIS A001383:
$$1, 1, 2, 4, 8, 15, 29, 53, 98, 177, 319, 565, \ldots$$
Finally $T_{\le 5}(z)$ yields the sequence
OEIS A001384:
$$1, 1, 2, 4, 9, 19, 42, 89, 191, 402, 847 \ldots$$
There are two interpretations here depending on whether the singleton
is supposed to have height zero or one. A relevant link is this
MSE link.
Remark. The complexity of the growth of the partition numbers
provides an idea of the difficulty of this problem, which is as
difficult if not more. Studying the OEIS entries in more detail it
does appear that simple recurrences for these can be computed.
Concerning the recurrence for the labeled case differentiation
produces
$$T'_{\le h}(z) = \exp T_{\le h-1}(z) +
z \exp T_{\le h-1}(z) \times T'_{\le h-1}(z)
\\ = \frac{1}{z} T_{\le h}(z) + T_{\le h}(z) \times T'_{\le h-1}(z).$$
Extracting coefficients here we get
$$n! [z^n] T'_{\le h}(z)
= T_{\le h, n+1} = n! [z^{n+1}] T_{\le h}(z)
+ n! [z^n ] T_{\le h}(z) \times T'_{\le h-1}(z)
\\ = \frac{1}{n+1} T_{\le h, n+1}
+ n! \sum_{q=1}^{n} [z^q] T_{\le h}(z) [z^{n-q}] T'_{\le h-1}(z)
\\ = \frac{1}{n+1} T_{\le h, n+1}
+ \sum_{q=1}^{n} {n\choose q} T_{\le h,q} T_{\le h-1, n-q+1}.$$
We obtain the closed form
$$\bbox[5px,border:2px solid #00A000]
{T_{\le h, n+1} = \frac{n+1}{n}
\sum_{q=1}^{n} {n\choose q} T_{\le h,q} T_{\le h-1, n-q+1}.}$$
The boundary condition here is $T_{\le h, 1} = 1$ and
$T_{0,1} = 0.$
Applying differentiation to the unlabeled case yields
$$T'_{\le h}(z) =
\exp\left(\sum_{l\ge 1} \frac{T_{\le h-1}(z^l)}{l}\right)
\\ + z\exp \left(\sum_{l\ge 1} \frac{T_{\le h-1}(z^l)}{l}\right)
\times \sum_{l\ge 1} T'_{\le h-1}(z^l) z^{l-1}
\\ = \frac{1}{z} T_{\le h}(z)
+ T_{\le h}(z) \times \sum_{l\ge 1} T'_{\le h-1}(z^l) z^{l-1}.$$
Extracting coefficients we find
$$[z^n] T'_{\le h}(z) = (n+1) T_{\le h, n+1}
\\ = [z^{n+1}] T_{\le h}(z)
+ \sum_{l=1}^n
\sum_{q=1}^n
[z^q] T_{\le h}(z) [z^{n-q}] z^{l-1} T'_{\le h-1}(z^l)
\\ = T_{\le h, n+1}
+ \sum_{l=1}^n
\sum_{q=1}^n
T_{\le h,q} [z^{n-q-l+1}] T'_{\le h-1}(z^l).$$
We obtain the closed form
$$T_{\le h, n+1} =
\frac{1}{n} \sum_{l=1}^n
\sum_{q=1}^n
T_{\le h,q} [[l|n-q+1]] ((n-q+1)/l) T_{\le h-1, (n-q+1)/l}$$
which simplifies to
$$\bbox[5px,border:2px solid #00A000]
{T_{\le h, n+1} =
\frac{1}{n}
\sum_{q=1}^n \sum_{l|n-q+1}
T_{\le h,q} ((n-q+1)/l) T_{\le h-1, (n-q+1)/l}.}$$
Boundary conditions are the same as in the labeled case.
This recurrence makes it possible to compute e.g. $T_{\le 6,n}$
which yields
$$1, 1, 2, 4, 9, 20, 47, 108, 252, 582,
\\ 1345, 3086, 7072, 16121,\ldots$$
which is OEIS A001385 and $T_{\le 7,n}$
which yields
$$1, 1, 2, 4, 9, 20, 48, 114, 278, 676,
\\ 1653, 4027, 9816, 23843,\ldots$$
which is OEIS A034823.
The Maple code to implement these two recurrences is as follows.
T :=
proc(h,m)
option remember;
local n;
if m=1 then if h=0 then return 0 else return 1 fi fi;
n := m-1;
(n+1)/n*
add(binomial(n,q)*T(h,q)*T(h-1,n-q+1), q=1..n);
end;
X :=
proc(h,m)
option remember;
local n;
if m=1 then if h=0 then return 0 else return 1 fi fi;
n := m-1;
1/n*
add(add(X(h,q)*
`if`(n-q+1 mod l =0, (n-q+1)/l*X(h-1,(n-q+1)/l), 0),
q=1..n), l=1..n);
end;
For Catalan trees let us recall the combinatorial class equation
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}
\mathcal{C} = \mathcal{Z} \times \textsc{SEQ}(\mathcal{C}).$$
which gives the functional equation
$$C(z) = z \frac{1}{1-C(z)}.$$
Marking the root degree we get
$$\mathcal{Z} \times
\mathcal{U}^0 \textsc{SEQ}_{=0}(\mathcal{C})
+ \mathcal{Z} \times
\mathcal{U}^1 \textsc{SEQ}_{=1}(\mathcal{C})
+ \mathcal{Z} \times
\mathcal{U}^2 \textsc{SEQ}_{=2}(\mathcal{C})
+ \cdots$$
This gives the generating function
$$R(z, u) = z \sum_{k\ge 0} u^k C(z)^k
= z \frac{1}{1-u C(z)}.$$
As a sanity check put $u=1$ to get $z/(1-C(z))$ which is correct.
Differentiate and put $u=1$ to obtain
$$\left. \frac{\partial}{\partial u} R(z, u) \right|_{u=1}
= \left. \frac{z}{(1-uC(z))^2} C(z) \right|_{u=1}
= \frac{1}{z} C(z)^2 C(z) = \frac{1}{z} C(z)^3.$$
Now using Lagrange inversion by residues (Egorychev method) we seek
$$[z^n] \frac{1}{z} C(z)^3 = [z^{n+1}] C(z)^3
= \frac{1}{n+1} [z^n] 3 C(z)^2 C'(z)
\\ = \frac{1}{n+1} \;\underset{z}{\mathrm{res}}\;
\frac{1}{z^{n+1}} 3 C(z)^2 C'(z).$$
Next put $C(z)=w$ to get with $z= w (1-w)$
$$\frac{1}{n+1} \;\underset{w}{\mathrm{res}}\;
\frac{1}{w^{n+1} (1-w)^{n+1}} 3w^2
= \frac{3}{n+1} [w^{n-2}] \frac{1}{(1-w)^{n+1}}
= \frac{3}{n+1} {n+n-2\choose n}
\\ = \frac{3}{n+1} {2n-2\choose n}.$$
It follows that the average degree of the root is
$$\frac{3}{n+1} {2n-2\choose n} n {2n-2\choose n-1}^{-1}
= \frac{3}{n+1} \frac{(n-1)!^2}{(n-1)! (n-2)!}$$
or
$$\bbox[5px,border:2px solid #00A000]{
3\frac{n-1}{n+1}.}$$
This result was verified with the combstruct package, see below.
with(combstruct);
with(combinat);
trees :=
proc(n)
option remember;
local spec, k, res;
res := 0;
for k to n-1 do
spec := { T1= Prod(Z, Sequence(T1)),
Z=Atom,
T2 = Prod(Z, Sequence(T1, card=k))};
res := res +
k*count([T2, spec, unlabeled], size=n);
od;
res/(1/n*binomial(2*n-2,n-1));
end;
Best Answer
You already showed that $[z^n]t(z) = \frac{1}{n}\binom{2(n-1)}{n-1}$ is the $(n-1)$st Catalan number. The number of labeled rooted plane trees is just this multiplied by $n!$ (accounting for the number of ways to label the vertices).
Using the symbolic method, this can be explained by noting that the GF for the labeled version also satisfies $\mathcal{T} = \mathcal{Z} \times \text{Seq}(\mathcal{T})$ (where the product and Sequence operations are for labeled GFs), so the GF is exactly the same as in the unlabeled case, so $n![z^n]t(z) = n! \frac{1}{n} \binom{2(n-1)}{n-1}$.