Number of all labeled, unordered rooted trees with $n$ vertices and $k$ leaves.

This answer explores the ideas by @vonbrand, using a form of Lagrange inversion to extract coefficients. I would write the species equation slightly differently, namely as

$$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{T} = \mathcal{Z} \times \mathcal{Y} + \mathcal{Z}\times \textsc{SET}_{\ge 1}(\mathcal{T}).$$

This translates into the functional equation $$T(z) = zy + z(\exp T(z) - 1) = z(-1+y+\exp T(z)).$$

We are interested in extracting coefficients as per the following equation: $$n! [z^n] T(z) = \frac{n!}{2\pi i}\int_{|z|=\epsilon} \frac{1}{z^{n+1}} T(z) dz.$$ Put $w=T(z)$ so that $$z = \frac{w}{-1+y+\exp w}$$ and $$dz = \left(\frac{1}{-1+y+\exp w} - \frac{w\exp w}{(-1+y+\exp w)^2}\right)\; dw.$$ This map fixes the origin, so that we may write $$\frac{n!}{2\pi i}\int_{|w|=\epsilon} \frac{(-1+y+\exp w)^{n+1}}{w^{n+1}} \times w \\ \times \left(\frac{1}{-1+y+\exp w} - \frac{w\exp w}{(-1+y+\exp w)^2}\right) dw.$$ We treat the two parenthesised terms in turn. The first one gives $$\frac{(-1+y+\exp w)^n}{w^n} = \frac{1}{w^n} \sum_{q=0}^n {n\choose q} (y-1)^{n-q} \exp(wq)$$ This produces the residue $$\sum_{q=0}^n {n\choose q} (y-1)^{n-q} \frac{q^{n-1}}{(n-1)!}.$$

The second one gives $$\exp w\frac{(-1+y+\exp w)^{n-1}}{w^{n-1}} = \frac{1}{w^{n-1}} \sum_{q=0}^{n-1} {n-1\choose q} (y-1)^{n-1-q} \exp(w(q+1))$$ This produces the residue $$\sum_{q=0}^{n-1} {n-1\choose q} (y-1)^{n-1-q} \frac{(q+1)^{n-2}}{(n-2)!}.$$

Let us check the case $y=1$ before we proceed to extract coefficients. This should give the total count of trees on $n$ nodes. Only the term for $q=n$ contributes from the first term, giving $$n! \times \frac{n^{n-1}}{(n-1)!} = n^n.$$ From the second term we get with $q=n-1$ $$n! \times \frac{n^{n-2}}{(n-2)!} = (n-1) \times n^{n-1}.$$ The difference of these two is $$n^{n-1}$$ which is the right value.

We return to the process of extracting coefficients from this and multiplying by $n!$ we get for the number of trees with $n$ nodes and $k$ leaves that $$n \times \sum_{q=0}^n {n\choose q} {n-q\choose k} (-1)^{n-q-k} q^{n-1} \\- n(n-1) \times \sum_{q=0}^{n-1} {n-1\choose q} {n-q-1\choose k} (-1)^{n-q-1-k} (q+1)^{n-2}.$$

Now with the following Maple definitions

c1 := proc(n,k)
   n*add(binomial(n,q)*binomial(n-q,k)*
         (-1)^(n-q-k)*q^(n-1), q=0..n);
end;

c2 := proc(n,k)
   n*(n-1)*add(binomial(n-1,q)*binomial(n-1-q,k)*
               (-1)^(n-1-q-k)*(q+1)^(n-2), q=0..n);
end;

c := (n, k) -> c1(n,k)-c2(n,k);

we get the following sequence for $k=1$ starting at $n=2:$ $$2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800,\ldots$$ which is correct and counts the number of rooted paths with a leaf at the bottom.

For $k=2$ starting at $n=3$ we get $$3, 36, 360, 3600, 37800, 423360, 5080320, 65318400, 898128000,\ldots$$ which points us to OEIS A055303 where we find that we have the right values.

Next try $k=3$ starting at $n=4$ to get $$4, 140, 3000, 54600, 940800, 16087680, 279417600, 4989600000, 92207808000,\ldots$$ which points us to OEIS A055304, once more confirming the correctness of the above derivation.

The case $k=4$ starting at $n=5$ yields $$5, 450, 18900, 588000, 15876000, 400075200, 9779616000, 237105792000, 5779453680000,\ldots$$ which points us to OEIS A055305. This too is confirmed.

The last one we will check is $k=5$ starting at $n=6$, which yields $$6, 1302, 101136, 5143824, 210198240, 7593173280, 255415628160, 8252203639680, 261173083691520,\ldots$$ which points to OEIS A055306, good news once again.

Additional OEIS consultation shows that even $k=10$ is listed there and there is a master table at OEIS A055302.

Addendum. The above formula for trees on $n$ nodes and with $k$ leaves admits additional simplification. Both terms are convolutions of exponential generating functions $A(z)$ and $B(z)$.

The first term has $$A(z) = \sum_{m\ge 0} {m\choose k} (-1)^{m-k} \frac{z^m}{m!} = \sum_{m\ge k} {m\choose k} (-1)^{m-k} \frac{z^m}{m!} \\= \frac{1}{k!} \sum_{m\ge k} \frac{(-1)^{m-k}}{(m-k)!} z^m = \frac{z^k}{k!} \exp(-z)$$ and $$B(z) = \sum_{m\ge 0} m^{n-1} \frac{z^m}{m!} = \exp(z) \sum_{q=1}^{n-1} {n-1\brace q} z^q$$ as is easily proved by induction.

This gives the following closed form for the first term: $$n \times n! [z^n] \frac{z^k}{k!} \sum_{q=1}^{n-1} {n-1\brace q} z^q = n \times n! [z^{n-k}] \frac{1}{k!} \sum_{q=1}^{n-1} {n-1\brace q} z^q \\= n \times \frac{n!}{k!} {n-1\brace n-k}.$$

The second term is closely related to the first, being evaluated at $n-1$ instead of $n.$ The function $A(z)$ is the same and $B(z)$ is $$B(z) = \sum_{m\ge 0} (m+1)^{n-2} \frac{z^m}{m!} = \exp(z) \sum_{q=1}^{n-1} {n-1\brace q} z^{q-1}$$ as is once more easily proved by induction.

This gives the following closed form for the second term: $$n(n-1) \times (n-1)! [z^{n-1}] \frac{z^k}{k!} \sum_{q=1}^{n-1} {n-1\brace q} z^{q-1} \\= n(n-1) \times (n-1)! [z^{n-1-k}] \frac{1}{k!} \sum_{q=1}^{n-1} {n-1\brace q} z^{q-1} \\ = n(n-1) \frac{(n-1)!}{k!} {n-1\brace n-k}.$$

The conclusion is that the number of trees on $n$ nodes and with $k$ leaves is given by $$\left(n - (n-1)\right) \frac{n!}{k!} {n-1\brace n-k} = \frac{n!}{k!} {n-1\brace n-k}.$$

Post scriptum. The simplest evaluation of the two $B(z)$ is in fact not by induction but uses a basic Stirling number identity. Suppose we seek $$\sum_{m\ge 0} m^{n-1} \frac{z^m}{m!}.$$ Now recall that $$m^{n-1} = \sum_{q=0}^{n-1} {n-1\brace q} m^{\underline q}.$$ Substitution yields $$\sum_{m\ge 0} \left(\sum_{q=0}^{n-1} {n-1\brace q} m^{\underline q}\right) \frac{z^m}{m!}$$ which is $$\sum_{q=0}^{n-1} {n-1\brace q} \sum_{m\ge 0} m^{\underline{q}} \frac{z^m}{m!} = \sum_{q=0}^{n-1} {n-1\brace q} \sum_{m\ge q} m^{\underline{q}} \frac{z^m}{m!} \\= \sum_{q=0}^{n-1} {n-1\brace q} \sum_{m\ge q} \frac{z^m}{(m-q)!} = \sum_{q=0}^{n-1} {n-1\brace q} z^q \exp(z)$$ and we have the result.

Use the analytic method. Your class is a root connected to a non-empty set of trees, or a leaf. Use $\mathcal{Z}$ (and $z$) for inner nodes, $\mathcal{Y}$ (and $y$) for leaves; use $\mathcal{E}$ for the class with one empty object: $$ \mathcal{T} = \mathcal{Z} \star (\mathfrak{S}(\mathcal{T}) \smallsetminus \mathcal{E}) + \mathcal{Z} \mathcal{Y} $$ This translates to: $$ T(z, y) = z (e^{T(z, y)} - 1) + z y $$ Just need to get $T(z, y)$ (or the coefficients) out of this...

I present an improvement over the first answer by way of enrichment. We start with the combinatorial class

$$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{T} = \mathcal{Z} \times \mathcal{Y} + \mathcal{Z}\times \textsc{SET}_{\ge 1}(\mathcal{T}).$$

This translates into the functional equation $$T(z) = zy + z\times(\exp T(z) - 1) = z\times(-1+y+\exp T(z))$$ or

$$z = \frac{T(z)}{-1+y+\exp T(z)}.$$

Writing

$$T(z) = \sum_{n\ge 1} T_n(y) \frac{z^n}{n!}$$

we are interested in extracting coefficients as per the Cauchy Coefficient Formula

$$\frac{1}{(n-1)!} T_n(y) = [z^{n-1}] T'(z) = \frac{1}{2\pi i}\int_{|z|=\epsilon} \frac{1}{z^{n}} T'(z) \; dz.$$

Now we put $w=T(z)$ so that $dw = T'(z) \; dz$ and use the functional equation to obtain

$$\frac{1}{2\pi i} \int_{|w|=\gamma} \frac{(-1+y+\exp w)^n}{w^{n}} \; dw.$$

Extracting the coefficient on $[y^k]$ we find

$$[y^k] \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{(-1+y+\exp w)^n}{w^{n}} \; dw \\ = \frac{1}{2\pi i} \int_{|w|=\gamma} {n\choose k} \frac{(\exp(w)-1)^{n-k}}{w^{n}} \; dw.$$

This yields

$$[y^k] T_n(y) = (n-1)! {n\choose k} [w^{n-1}] (\exp(w)-1)^{n-k} \\ = [y^k] T_n(y) = \frac{n!}{k!} (n-1)! [w^{n-1}] \frac{(\exp(w)-1)^{n-k}}{(n-k)!}.$$

We recognize the EGF of the Stirling numbers of the second kind, getting

$$\bbox[5px,border:2px solid #00A000]{ \frac{n!}{k!} {n-1\brace n-k}.}$$

This is from Flajolet & Sedgewick, Analytic Combinatorics, page 732, Lagrange Inversion.

As a sanity check we find for the number of all rooted trees

$$\sum_{k=1}^n (n-1)! {n\choose k} [w^{n-1}] (\exp(w)-1)^{n-k} \\ = (n-1)! [w^{n-1}] \sum_{k=1}^n {n\choose k} (\exp(w)-1)^{n-k}.$$

Now for $k=0$ we get

$$(n-1)! [w^{n-1}] (\exp(w)-1)^n = 0$$

because $\exp(w)-1 = w + \cdots.$ Hence we may continue with

$$(n-1)! [w^{n-1}] \sum_{k=0}^n {n\choose k} (\exp(w)-1)^{n-k} = (n-1)! [w^{n-1}] \exp(nw) \\ = (n-1)! \frac{n^{n-1}}{(n-1)!} = n^{n-1}$$

and the check goes through (Cayley's formula).