I am sending some preliminary, as yet not fine-tuned and optimized
obervations to generate activity on this question. I am sure these
admit improvements (the time and space complexity of the algorithm is
poor).
If I have decoded the problem statement correctly we are enumerating
labeled endofunctions with no fixed points as appeared at this MSE
link so that these
pseudo-forests are directed.
The species $\mathcal{Q}$ under consideration is
$$\mathcal{Q}
= \mathfrak{P}(\mathfrak{C}_{=2}(\mathcal{T})
+ \mathfrak{C}_{=3}(\mathcal{T})
+ \mathfrak{C}_{=4}(\mathcal{T}) + \cdots).$$
where $\mathcal{T}$ represents labeled rooted trees with EGF $T(z)$
and functional equation $T(z) = z\exp T(z),$ the labeled tree
function.
It follows that the distribution of the maxima of the tree sizes
(where unicyclic components are counted by the number of nodes) on $n$
nodes can be found by considering the species
$$\mathcal{Q}_{\le n}
= \mathfrak{P}(\mathfrak{C}_{=2}(\mathcal{T}_{\le n})
+ \mathfrak{C}_{=3}(\mathcal{T}_{\le n})
+ \mathfrak{C}_{=4}(\mathcal{T}_{\le n}) + \cdots
+ \mathfrak{C}_{=n}(\mathcal{T}_{\le n}))$$
and marking the cycles with a variable $\mathcal{V}_q$ indicating the
number of nodes in the component using the generating function
$$Q(z) = \exp
\left(\frac{T_{\le n}(z)^2}{2}
+ \frac{T_{\le n}(z)^3}{3}
+ \frac{T_{\le n}(z)^4}{4} + \cdots
+ \frac{T_{\le n}(z)^n}{n}\right).$$
This will produce the following distributions:
$${u}^{2},\\
8\,{u}^{3},\\
78\,{u}^{4}+3\,{u}^{2},\\
944\,{u}^{5}+80\,{u}^{3},\\
13800\,{u}^{6}+1170\,{u}^{4}+640\,{u}^{3}+15\,{u}^{2},\\
237432\,{u}^{7}+19824\,{u}^{5}+21840\,{u}^{4}+840\,{u}^{3},\\
4708144\,{u}^{8}+386400\,{u}^{6}+422912\,{u}^{5}+229320\,{u}^{4}
\\+17920\,{u}^{3}+105\,{u}^{2},\\
105822432\,{u}^{9}+8547552\,{u}^{7}+9273600\,{u}^{6}+9634464\,{u}^{5}\\
+786240\,{u}^{4}+153440\,{u}^{3},\ldots$$
This gives the following for the expected maximum tree size:
$$2., 3., 3.925925926, 4.843750000, 5.723520000, 6.612311385,
\\ 7.471092584, 8.342072010, 9.189007167, 10.04727275,
\\ 10.88589292, 11.73525388, 12.56739638, 13.40959924,\ldots$$
which would seem to indicate that most of these consist of one
connected component.
The following Maple code was used to compute these values. It took
$58$ seconds to compute the distribution for $n=35,$
which represents $$399725722782532944388077717044552088857010024925364224$$
pseudoforests ($34^{35}$).
gf_le :=
proc(n)
option remember;
local Tle, term, res, p;
Tle := add(q^(q-1)*z^q/q!, q=1..n);
res := 0;
for term in expand(add(Tle^q/q, q=2..n)) do
p := degree(term, z);
res := res + v[p]*coeff(term, z, p)*z^p;
od;
exp(res);
end;
gf :=
proc(n)
option remember;
local res, cft, vs, term;
res := 0; cft := n!*coeftayl(gf_le(n), z=0, n);
if not type(cft, `+`) then
cft := [cft];
fi;
for term in cft do
vs := indets(term);
res := res +
lcoeff(term)*u^max(map(t->op(1, t), vs));
od;
res;
end;