Here are some computational aspects of this problem that can serve as
a basis for further exploration. We present relevant generating
functions and recurrence relations for unlabelled trees and
endofunctions. These basic data of course cannot take the place of a
proper investigation using combinatorial methods.
We have by inspection that the combinatorial class $\mathcal{M}$ under
consideration here is (notation from Analytic Combinatorics)
$$\mathcal{M} = \def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}\textsc{MSET}(\textsc{CYC}(\mathcal{T}))$$
with $\mathcal{T}$ being the combinatorial class of rooted trees.
Therefore the study of $\mathcal{M}$ requires the study of
$\mathcal{T}.$
We have for rooted trees that
$$\mathcal{T} = \mathcal{Z} \times \textsc{MSET}(\mathcal{T}).$$
This translates to the functional equation
$$T(z) = z \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right).$$
Differentiate to get a recurrence for the number $T_n$ of rooted trees:
$$T'(z) = \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right)
+ z \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right)
\left(\sum_{l\ge 1} \frac{T'(z^l)}{l} \times l z^{l-1} \right)
\\ = \frac{T(z)}{z}
+ T(z) \left(\sum_{l\ge 1} T'(z^l) z^{l-1} \right) .$$
so that
$$z T'(z) = T(z) + T(z) \left(\sum_{l\ge 1} T'(z^l) z^l \right)
= T(z) + T(z)
\left(\sum_{l\ge 1} z^l \sum_{m\ge 1} m T_m z^{(m-1)l}\right)$$
which finally yields
$$z T'(z) = T(z) + T(z)
\left(\sum_{l\ge 1} \sum_{m\ge 1} m T_m z^{ml}\right).$$
Extracting coefficients from this we get
$$n T_n = T_n + [z^n] T(z)
\left(\sum_{l\ge 1} \sum_{m\ge 1} m T_m z^{ml}\right)
= T_n + \sum_{l=1}^{n-1}
\sum_{m=1}^{\lfloor(n-1)/l\rfloor} m T_m T_{n-ml}$$
so that
$$T_n = \frac{1}{n-1}
\sum_{l=1}^{n-1}
\sum_{m=1}^{\lfloor(n-1)/l\rfloor} m T_m T_{n-ml}$$
with $T_0 = 0$ and $T_1 = 1.$
This gives the sequence
$$1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766,
12486, 32973, 87811, \ldots$$
which is OEIS A000081
where the above derivation is confirmed.
Now return to the combinatorial class $\mathcal{Q}$ and compute the generating
function of the operator $\textsc{MSET}(\textsc{CYC}(\cdot)).$
Using the cycle index of the cyclic group we easily derive that the
generating function of the cycle operator being applied to some
combinatorial class $\mathcal{A}$ is
$$\sum_{k\ge 1} \frac{\varphi(k)}{k}
\log\frac{1}{1-A(z^k)}.$$
Hence the generating function for
$\textsc{MSET}(\textsc{CYC}(\cdot))$ is
$$\exp\left(
\sum_{l\ge 1}\frac{1}{l} \sum_{k\ge 1} \frac{\varphi(k)}{k}
\log\frac{1}{1-A(z^{kl})}\right).$$
This is
$$\exp
\left(\sum_{m\ge 1} \log\frac{1}{1-A(z^m)}
\sum_{k|m} \frac{\varphi(k)}{m}\right)
\\ = \exp
\left(\sum_{m\ge 1} \log\frac{1}{1-A(z^m)}\right)
= \prod_{m\ge 1} \frac{1}{1-A(z^m)}.$$
It follows that the generating function $Q(z)$ of $\mathcal{Q}$ in
terms of $T(z)$ is given by
$$Q(z) = \prod_{m\ge 1} \frac{1}{1-T(z^m)}.$$
Differentiate to obtain a recurrence from this:
$$Q'(z) = \prod_{m\ge 1} \frac{1}{1-T(z^m)}
\times \sum_{m\ge 1} (1-T(z^m))
\frac{T'(z^m) \times m z^{m-1}}{(1-T(z^m))^2}
\\ = Q(z) \times \sum_{m\ge 1}
\frac{T'(z^m) \times m z^{m-1}}{1-T(z^m)}.$$
Extracting coefficients from this we obtain
$$(n+1) Q_{n+1} =
\sum_{q=0}^n Q_{n-q}
[z^q] \sum_{m\ge 1} \frac{T'(z^m) \times m z^{m-1}}{1-T(z^m)}
\\ = \sum_{q=0}^n Q_{n-q}
[z^q] \sum_{m\ge 1} \frac{m}{z}
\frac{\sum_{p\ge 1} p T_p z^{m (p-1)} \times z^m}{1-T(z^m)}
\\ = \sum_{q=0}^n Q_{n-q}
\sum_{m=1}^{q+1} m [z^{q+1}]
\frac{\sum_{p\ge 1} p T_p z^{mp}}{1-T(z^m)}.$$
Now here we need to digress for a moment to find a recurrence for the
coefficients $P_n$ of
$$P(z) = \frac{1}{1-T(z)}.$$
Differentiating yields
$$P'(z) = \frac{T'(z)}{(1-T(z))^2} = P(z)^2 T'(z).$$
Extracting coefficients we obtain
$$(n+1) P_{n+1}
= \sum_{q=0}^n \left(\sum_{p=0}^{n-q} P_p P_{n-q-p}\right)
[z^q] T'(z)
\\ = \sum_{q=0}^n (q+1) T_{q+1}
\left(\sum_{p=0}^{n-q} P_p P_{n-q-p}\right)$$
so that
$$P_n =
\frac{1}{n}
\sum_{q=0}^{n-1} (q+1) T_{q+1}
\left(\sum_{p=0}^{n-1-q} P_p P_{n-1-q-p}\right)$$
with $P_0 = 1.$
This yields the sequence
$$1, 2, 5, 13, 35, 95, 262, 727, 2033, 5714, 16136,
45733, 130046, 370803,\ldots$$
which is OEIS A000107 where a better
recurrence can also be found.
Returning to $Q_{n+1}$ and using an Iverson bracket we obtain
$$\sum_{q=0}^n Q_{n-q}
\sum_{m=1}^{q+1} m [[m|q+1]]
\sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}.$$
which yields the formula
$$Q_n = \frac{1}{n}
\sum_{q=0}^{n-1} Q_{n-1-q}
\sum_{m=1}^{q+1} m [[m|q+1]]
\sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}$$
where $Q_0 = Q_1 = 1.$
An alternate form of this is
$$Q_n = \frac{1}{n}
\sum_{q=0}^{n-1} Q_{n-1-q}
\sum_{m|q+1}
m \sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}$$
again where $Q_0 = Q_1 = 1.$
This finally gives the sequence
$$1, 3, 7, 19, 47, 130, 343, 951, 2615, 7318, 20491, 57903, 163898,
\\ 466199, 1328993,\ldots$$
which points us to OEIS A001372 where additional material awaits and the results of the above calculation are validated.
The following Maple code implements these three sequences.
with(numtheory, divisors);
T :=
proc(n)
option remember;
if n<=1 then return n fi;
1/(n-1)*add(add(m*T(m)*T(n-m*l),
m=1..floor((n-1)/l)), l=1..n-1);
end;
P :=
proc(n)
option remember;
if n=0 then return 1 fi;
1/n*add((q+1)*T(q+1)*
add(P(p)*P(n-1-q-p), p=0..n-1-q), q=0..n-1);
end;
Q :=
proc(n)
option remember;
local res, iv, q, m;
if n<=1 then return 1 fi;
1/n*add(Q(n-1-q)*
add(m*add(p*T(p)*P((q+1)/m-p), p=1..(q+1)/m),
m in divisors(q+1)), q=0..n-1);
end;
Addendum.
Note that just for curiosity's sake we can do the labelled
case es well, which is very simple. Labelled trees are given by
$$\mathcal{T} =
\mathcal{Z} \times \textsc{SET}(\mathcal{T})$$
which gives the functional equation
$$T(z) = z \exp T(z).$$
The combinatorial class of endofunctions is then given by
$$\mathcal{Q} = \textsc{SET}(\textsc{CYC}(\mathcal{T})).$$
Translating to generating functions we get
$$Q(z) = \exp \log \frac{1}{1-T(z)} =
\frac{1}{1-T(z)}.$$
Extracting coefficients via Lagrange inversion we have
$$Q_n
= n! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{1}{1-T(z)} dz.$$
Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and
$dz = \exp(-w) - w\exp(-w)$
to get
$$n! \frac{1}{2\pi i}
\int_{|w|=\epsilon}
\frac{\exp(w(n+1))}{w^{n+1}}
\frac{1}{1-w} (\exp(-w) - w\exp(-w)) dw
\\ = n! \frac{1}{2\pi i}
\int_{|w|=\epsilon}
\frac{\exp(wn)}{w^{n+1}}
\frac{1}{1-w} (1 - w) dw
\\ = n! \frac{1}{2\pi i}
\int_{|w|=\epsilon}
\frac{\exp(wn)}{w^{n+1}} dw.$$
But we have
$$n! [w^n] \exp(wn) = n! \times \frac{n^n}{n!} = n^n$$
and we have just verified the count of the number of
endofunctions from basic algebra.
All of the above (labelled and unlabelled) is already present in some form or other in Harary and Palmer's Graphical Enumeration.
Remark I. Care should be taken not to confuse the labelled and unlabelled tree functions.
Remark II. Care should furthermore be taken in the proper usage of sets and multisets in the labelled and the unlabelled case. We distinguish multisets and sets in the unlabelled case but in the labeled case the cartesian product that underlies the construction of all combinatorial operators involves re-labeling that uniquely identifies the two components by the choice of labels and there can never be multisets (operator $\textsc{MSET}$), just sets (operator $\textsc{SET}$).