(PARI) \\ See links in
A339645 for combinatorial species functions.
edges(v) = {sum(i=2, #v, sum(j=1, i-1, gcd(v[i], v[j]))) + sum(i=1, #v, v[i]\2)}
graphsCycleIndex(n)={my(s=0); forpart(p=n, s+=permcount(p) * 2^edges(p) * sMonomial(p)); s/n!}
graphsSeries(n)={sum(k=0, n, graphsCycleIndex(k)*x^k) + O(x*x^n)}
cycleIndexSeries(n)={my(g=graphsSeries(n), gcr=sPoint(g)/g); x*sSolve( sLog( gcr/(x*sv(1)) ), gcr )}
{ my(N=15); Vec(OgfSeries(cycleIndexSeries(N)), -N) } \\
Andrew Howroyd, Dec 25 2020