Tanglegrams and functional equations of M. Somos Recent references on the matter at hand include, a lecture slide The Konvalinka-Amdeberhan conjecture
and plethystic inverses and a preprint on Counting tanglegrams with species by I. Gessel; the initial work On the enumeration of tanglegrams and tangled chains by S. Billey et al.
For each prime $p$, let's consider the family of (integral) sequences
$$a_p(n):=\sum_{\lambda\vdash n}\frac{1}{z_\lambda}\prod_{i=2}^{\ell(\lambda)}(p\lambda_i+\cdots+p\lambda_{\ell(\lambda)}+1);$$
where the sum runs through all $p$-ary partitions $\lambda$, of $n$, and $\ell(\lambda)$ is the length of the partition and the numbers $z_{\lambda}$ are well-known since the number of permutations in $\mathfrak{S}_n$ with cycle type $\lambda$ is computed by $\frac{n!}{z_{\lambda}}$.
On the other hand, Michael Somos proposed the functional equations $$
x\cdot A_p(x)^p=\frac{x\cdot A_p(x^p)}{1-p\cdot x\cdot A_p(x^p)}
$$ and a couple of these are listed on OEIS as A085748 and A091190. However, there are not enough interpretations attached to $A_p(x)$ or its coefficients on OEIS.
After some experimentation, I wish to ask:

QUESTION. Can the following be justified or refuted?
$$A_p(x)=\sum_{n\geq0}a_p(n)x^n.$$

Remark. Observe that neither $A_p(x)$ nor $a_p(n)$ have been found to output integers, unless $p$ is a prime.
 A: In my lecture slides that Tewodros cites, I studied symmetric functions that I called $g_m$, where $m$ is a positive integer, that have constant term 1 and satisfy
$$ -L_m[g_m] = p_1,$$
where $L_m$ is the Lyndon symmetric function given by
$$L_m = \frac{1}{m} \sum_{d\mid m} \mu(d) p_d^{m/d}.$$
Here $\mu$ is the Möbius function, the $p_i$ are power sum symmetric functions, and $L_m[g_m]$ is the plethysm of $L_m$ and $g_m$. I showed that $g_m$ is an integral symmetric function; i.e., its coefficients in the underlying variables are integers, and that if $m$ is a power of the prime $q$ then for any $\alpha$,
\begin{multline*}
\quad
g_m^{-\alpha} = 1+
\sum_{n=1}^\infty \sum_{\lambda}\frac{p_{\lambda}}{z_\lambda}
\times
\alpha\prod_{j=2}^{l(\lambda)} (m\lambda_j+m\lambda_{j+1}+\cdots +m\lambda_{l(\lambda)}+\alpha),\quad\tag{1}
\end{multline*}
where the sum on $\lambda$ is over all partitions of $n$ in which every part is a power of $q$. (I'm using $q$ here instead of $p$ to avoid confusion with the power sum symmetric functions.) This shows that if $\alpha$ is positive and $m$ is a power of a prime then all coefficients of $g_m^{-\alpha}$ are positive. I don't know that this is true if $m$ is not a power of a prime, though I suspect that it is.
Now let us specialize the symmetric functions (in the variables $x_1, x_2,\dots$) by setting $x_1=x$ and $x_i=0$ for $i>0$, or equivalently, $p_i =x^i$ for all $i$. Let $G_m(x)$ be the image of $g_m$ under this specialization. Then $G_m(x)$ satisfies the functional equation
$$-\frac{1}{m}\sum_{d\mid m}\mu(d) G_m(x^d)^{m/d} = x.$$
If $m$ is a power of the prime $p$, this may be written
$$G_m(x^p)^{m/p} -G_m(x)^m = mx$$
and it's easy to check that if $m=p$ then $A_p(x) =1/G_p(x)$, where $A_p(x)$ is as in the original question. Moreover, it follows from $(1)$ that if $m$ is a power of the prime $p$ then for any $\alpha$ we have
\begin{multline*}
\quad
G_m(x)^{-\alpha} = 1+
\sum_{n=1}^\infty x^n
\sum_{\lambda}\frac{\alpha}{z_\lambda}\, 
\times\prod_{j=2}^{l(\lambda)} (m\lambda_j+m\lambda_{j+1}+\cdots +m\lambda_{l(\lambda)}+\alpha),\quad\tag{2}
\end{multline*}
where the sum on $\lambda$ is over all $p$-ary partitions of $n$. The case $m=p$, $\alpha=1$ answers the OP's questions in the affirmative.
In all cases $G_m(x)$ has integer coefficients, but I can only prove that $G_m(x)^{-1}$ has positive coefficients when $m$ is a prime power, when it follows from $(2)$, though this is likely true in all cases. The stronger statement that $1-G_m(x)$ has positive coefficients also follows from $(2)$ when $m$ is prime power.
I am working on a paper with detailed proofs of these formulas, but it will take a while.
