Follow-up on Rupinski's and Chapoton's observations:
To nail down the identification of the $p_n(x)$ with the cycle index polynomials for $S_n$ (or the partition polynomials of the refined Stirling numbers of the first kind A036039), look at the Taylor series rep of the digamma operator for the raising / creation operator for the $p_n(z)$ basis
$$R_z = z-\Psi(1+D_z) = z+\gamma + \sum_{n=1}^{\infty } (-1)^n\zeta (n+1)D_z^n.$$
This is precisely the raising operator for the cycle index polynomials as presented on page 23 of Lagrange à la Lah Part I with $c_1=z+\gamma=p_1(x)$ and $c_{n+1}=(-1)^n\zeta(n+1)$ for $n>0$
$$D^{-1}_{c_1}= :\frac{c_{.}}{1-c_{.}D_{c_1}}: = c_1+\sum_{n=1}^{\infty } c_{n+1}D_{c_1}^n.$$
Alternatively, the Newton identities extrapolated to an entire function as an infinite order polynomial using the Weierstrass factorization maneuver can be applied to see the connections to the power and elementary symmetric polynomial formalism:
$$\exp\left (-\beta p_{.}(z)\right )=\frac{\exp\left (-\beta z \right )}{\left (-\beta \right )!}=\exp\left (-\beta(z+\gamma) \right )\prod_{k=1}^{\infty }\left ( 1-\frac{\beta}{k} \right )\exp\left (\frac{\beta}{k} \right )$$
$$=\exp\left [-(z+\gamma)\beta -\sum_{k=2}^{\infty } \frac{\zeta (k)\beta ^k}{k} \right ]=\exp\left [ :ln(1-a\beta ) :\right ]$$ where $a^1=a_{1}=(z+\gamma)$ and $a^k=a_k=\zeta(k)$ for $k>1$, but this is precisely the umbral form of the e.g.f. for the cycle index polynomials (mod signs).
(Also there are connections to rational zeta series.)
Update (Nov. 16, 2012): The generating series appears on pg. 58 in "Hodge theoretic aspects of mirror symmetry" by L. Katzarkov, M. Kontsevich, and T. Pantev (following Lu's references).
