As I started working on this using Burnside's lemma, I might as well show what I've got. It is easy to make a mistake here.

Marko Riedel cast a serious doubt on the result. Lo and behold! It turned out that I had made an error in calculating that fraction. Because the erroneous sum still produced an integer I didn't catch it myself.

As described in loc.cit the 24 rotations of the cube split into 5 conjugacy classes. For each class we need to know what kind of orbits a representative in that class (or rather, the cyclic group generated by the representative) has among the edges. We also need to know the cardinality of the class.

$A$: The class of the identity mapping. Its action on the edges has 12 singleton orbits, and $|A|=1$.

$B$: The class of 90 degree rotations. Its action partitions the edges into 3 orbits with 4 edges in each orbit. $|B|=6.$

$C$: The class of 180 degree rotations about an axis connecting the midpoints of an opposite pair of faces of the cube (the squares of elements of class $B$). The edges are partitioned into 6 pairs. $|C|=3.$

$D$: The class of 120 degree rotations about an axis connecting two opposite vertices. The edges are partitioned to 4 orbits, each with 3 edges. $|D|=8.$

$E$: The class of 180 degree rotations about an axis connecting the midpoints of two parallel but opposite edges. The edges are partitioned into 7 orbits: 2 singletons (the ones intersecting the axis) and 5 pairs. $|E|=6.$

The group of color permutations is isomorphic to $S_3$ and has three conjugacy classes: $a$= the class of the identity mapping, $b$= the class of three 2-cycles (this time swapping two colors and leaving the third color as it is), $c$= the class of two 3-cycles.

The cartesian product of order $24\cdot6=144$ then acts on the set $X$ of colorings. The following table lists the number of fixed points in $X$ that an element of the listed conjugacy class has.

$$ \begin{array}{c|c|ccc|c|c} \text{class}&\text{size}&|X^g|&&\text{class}&\text{size}&|X^g|\\ \hline Aa&1&3^{12}&&Cc&6&0\\ Ab&3&1&&Da&8&3^4\\ Ac&2&0&&Db&24&1\\ Ba&6&3^3&&Dc&16&3^4\\ Bb&18&3^3&&Ea&6&3^7\\ Bc&12&0&&Eb&18&3^5\\ Ca&3&3^6&&Ec&12&0\\ Cb&9&3^6&&&& \end{array} $$ Feeding these numbers into the formula of Burnside's lemma tells that the number of orbits (=the number of distinguishable colorings under the prescribed rules) is $$ |X/G|=\frac{1\cdot3^{12}+3\cdot1+6\cdot3^3+18\cdot3^3+3\cdot3^6+9\cdot3^6+8\cdot3^4+24\cdot1+16\cdot3^4+6\cdot3^7+18\cdot3^5}{144}=3891. $$

How did I come up with these numbers. The details vary a bit case-by-case. As somewhat typical examples consider:

The class $Ab$: Here the edges stay at the same place, but two of the colors are interchanged. A coloring is a fixed point of this operation only if it is monochromatic (of the non-interchanged color). Thus the elements in this class have exactly one fixed point.

The class $Ba$: Here we want to calculate the number of colorings that remain unchanged, when we rotate the cube by 90 degrees, and don't do any repainting. For this to happen it is necessary and sufficient that any of the three groups of four edges permuted among themselves by the rotation must share the same color. We have the liberty of choosing the colors for all the orbits independently. Therefore there are $3^3$ such colorings.

The class $Bb$: Here the 3 groups of 4 edges undergo some repainting in addition to the 90 degree rotation. If the three colors are R,Y,G, and we are to repaint every yellow edge green, and every green edge yellow, then for each orbit the possible colorings are RRRR, YGYG and GYGY. In the latter two cases repainting cancels the effect of the rotation. So we have 3 choice for the three groups of 4 edges, and a total of $3^3$ colorings fixed under a group action of this type.

The class $Bc$: Here we can assume a color replacement scheme $Y\mapsto R\mapsto G\mapsto Y$. This permutation of order 3 just doesn't mesh with the rotation of order 4, so no colorings are fixed points for this operation.

The class $Dc$: Let's study the same color replacement scheme as with class $Bc$. The 12 edges are partitioned into 4 groups of 3 edges. The coloring of such a group has to be either $RGY$, $GYR$ or $YRG$ for the repainting to cancel the effect of the rotation. But the choices for the four orbits are independent from each other, so we get a total of $3^4$ colorings that are fixed under an operation of this type.

The class $Eb$: Assume color replacement scheme $R\mapsto R$, $Y\mapsto G\mapsto Y$. For a coloring to be a fixed point, the two edges that are rotated in place, must both be red. But for the 5 pairs of edgese that interchange their locations, the possible colorings are $RR$, $YG$ and $GY$ for the same reasons as with class $Bb$.

If you have trouble filling in the rest of the table, @ping me!


I would like to contribute some enrichment material to facilitate additional exploration of this problem. The most important observation is that what we have here is an instance of Power Group Enumeration, with the group acting on the slots where the $n$ colors are placed being the edge permutation group $E$ of the cube and the group acting on the colors being the symmetric group on $n$ elements $S_n$.

We can compute the number of configurations by Burnside's lemma which says to average the number of assignments fixed by the elements of the power group, which has $n!\times |E|!$ elements and $|E|=24$. But this number is easy to compute. Suppose we have a permutation $\alpha$ from $E$ and a permutation $\beta$ from $S_n.$ If we place the appropriate number of complete, directed and consecutive copies of a cycle from $\beta$ on a cycle from $\alpha$ then this assignment is fixed under the power group action, and this is possible iff the length of $\beta$ divides the length of $\alpha$ and there are as many assignments as the length of $\beta.$

To do this we first need the cycle index of $E$, which we now compute. There is the identity, which contributes $$a_1^{12}.$$ Rotations about an axis passing through a pairs of opposite vertices contribute $$4\times 2 a_3 ^4.$$ Rotations about an axis passing through opposite faces contribute $$3 \times (2 a_4^3 + a_2^6).$$ Finally rotations about an axis passing through midpoints of opposite edges contribute $$6\times a_1^2 a_2^5.$$

This gives the cycle index $$Z(E) = \frac{1}{24} \left(a_1^{12} + 8 a_3^4 + 6 a_4^3 + 3 a_2^6 + 6 a_1^2 a_2^5\right).$$

Let's verify this cycle index before we proceed. For edge colorings with $n$ colors where the colors are not being permuted we obtain the formula $$\frac{1}{24} \left(n^{12} + 8 n^4 + 6 n^3 + 3 n^6 + 6 n^7\right).$$ This gives the sequence $$1, 218, 22815, 703760, 10194250, 90775566, 576941778, 2863870080,\ldots$$ which points us to OEIS A060530, where we find that indeed we have the right cycle index.

Now the Burnside computation is best done with a CAS, here is the Maple code.

pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

cube_edge_cind :=
1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5);

cube_edge_colorings :=
proc(n)
option remember;
local idx_cols, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if n = 1 then return 1 fi;
    idx_cols := pet_cycleind_symm(n);

    res := 0;

    for term_a in cube_edge_cind do
        for term_b in idx_cols do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b*inst_b;
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;

This will produce the following sequence of values for colorings with $n$ colors with the symmetric group acting on the colors: $$1, 114, 3891, 29854, 87981, 143797, 170335, 177160,\ldots$$

Addendum Jan 29 2019. The above computation is easily adapted to produce a generating function, making it possible not only to count these colorings but also to classify them by the numbers of interchangeable colors used. The technique is the same (cover the cycles of a term from $Z(E)$ by cycles obtained from a permutation belonging to $S_n$), only this time we record the colors being used. With the colors being interchangeable we collect the contributions corresponding to the same partition of the twelve edges according to the colors into one term. This is the result for three colors:

$${P_{{1}}}^{12}+29\,{P_{{1}}}^{6}{P_{{2}}}^{6} +38\,{P_{{1}}}^{5}{P_{{2}}}^{7}+27\,{P_{{1}}}^{4}{P_{{2}}}^{8} +282\,{P_{{1}}}^{4}{P_{{2}}}^{4}{P_{{3}}}^{4} +13\,{P_{{1}}}^{3}{P_{{2}}}^{9} \\+1170\,{P_{{1}}}^{3}{P_{{2}}}^{4}{P_{{3}}}^{5} +412\,{P_{{1}}}^{3}{P_{{2}}}^{3}{P_{{3}}}^{6} +5\,{P_{{1}}}^{2}{P_{{2}}}^{10} +370\,{P_{{1}}}^{2}{P_{{2}}}^{5}{P_{{3}}}^{5} \\+600\,{P_{{1}}}^{2}{P_{{2}}}^{4}{P_{{3}}}^{6} +340\,{P_{{1}}}^{2}{P_{{2}}}^{3}{P_{{3}}}^{7} +77\,{P_{{1}}}^{2}{P_{{2}}}^{2}{P_{{3}}}^{8} +P_{{1}}{P_{{2}}}^{11} +236\,P_{{1}}{P_{{2}}}^{5}{P_{{3}}}^{6} \\+170\,P_{{1}}{P_{{2}}}^{4}{P_{{3}}}^{7} +85\,P_{{1}}{P_{{2}}}^{3}{P_{{3}}}^{8} +30\,P_{{1}}{P_{{2}}}^{2}{P_{{3}}}^{9} +5\,P_{{1}}P_{{2}}{P_{{3}}}^{10}.$$

The Maple code goes as follows.

with(combinat);

cube_edge_cind :=
1/24*(a[1]^12 + 8*a[3]^4 + 6*a[4]^3 + 3*a[2]^6 + 6*a[1]^2*a[2]^5);

pet_cycleind_symm :=
proc(n)
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;


pet_indets2rep :=
proc(ip)
local rep, var, deg, pos;

    rep := []; pos := 1;

    for var in indets(ip) do
        for deg to degree(ip, var) do
            rep :=
            [op(rep), [seq(s, s=pos..pos+op(1, var)-1)]];
            pos := pos + op(1, var);
        od;
    od;

    rep;
end;

cube_edge_colorings_gf :=
proc(n)
option remember;
local idx_cols, rep, res, term_a, term_b,
    v_a, v_b, inst_a, len_a, len_b, p, q,
    parts, alldeg, cols;

    if n = 1 then return P[1]^12 fi;
    idx_cols := pet_cycleind_symm(n);

    res := 0;

    for term_a in cube_edge_cind do
        for term_b in idx_cols do
            rep := pet_indets2rep(term_b);

            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in rep do
                    len_b := nops(v_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b*
                        mul(P[col], col in v_b)
                        ^(len_a/len_b);
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;

    parts := 0;

    for cols in expand(res) do
        alldeg :=
        sort(map(v -> degree(cols, v),
                 [op(indets(cols))]));

        parts := parts +
        lcoeff(cols)*
        mul(P[v]^alldeg[v], v=1..nops(alldeg));
    od;

    parts;
end;

cube_edge_colorings :=
proc(n)
option remember;
local gf, sl;

    gf := cube_edge_colorings_gf(n);
    sl := [seq(v = 1, v in indets(gf))];

    subs(sl, gf);
end;

I wrote up a how-to with several detailed examples, including counting the number of ways of coloring the faces of the cube, on my blog. That should get you started. I may be able to come back to this answer later to discuss edge colorings, but it's not very different.

Addendum: Jyrki has done a very thorough discussion of the fixed points of edge colorings under various actions, and I have nothing much to add to his excellent explanantion except that it is often only little more work to do the analysis in general for $N$ colors, and to get a polynomial in $N$, than it is to do it for some particular $N$, as Jyrki did.