If you take just one symbol with rotational symmetry of order 2 then there are 224 distinct cubes which can be made. A possible approach would be to assign symbol orientation and tab/slot at random to each of the 144 internal face pairs of the assembled $4\times 4\times 4$ puzzle; verify that the fully internal cubes are distinct; use bipartite matching to assign symbols to the external faces such that you get distinct cubes and no two corner cubes are interchangeable, and then try to verify the absence of phantom solutions with a SAT solver.

"I think (tho I haven't sat down to prove) that rotations will never fix a spanning tree." Counterexample when $d=5$: edge set $\{(-x,+z), (-y,+x), (-z,+y), (-v,-x), (-v,-y), (-v,-z), (-v,-w), (-w,+v), (+v,+w)\}$ is fixed by rotations which permute $(x, y, z)$.

I brute forced all spanning trees up to $d=5$, grouped them by the automorphisms of the hyperoctahedral graph, and get (group size, frequency) data $[(24, 6), (48, 5)]$ ($d=3$); $[(48, 4), (96, 4), (192, 77), (384, 176)]$ ($d=4$); $[(80, 2), (240, 6), (320, 2), (480, 60), (640, 2), (960, 147), (1280, 2), (1920, 1971), (3840, 7502)]$ ($d=5$). That the higher frequencies occur for group sizes $2^{-j}$ times the size of the automorphism group suggests that a careful account of trees fixed by $j$ reflections might get a reasonably tight estimate.

@BrendanMcKay, my bad: for "hypercube graph" substitute "hyperoctahedral graph".

From the main result of DeSplinter et. al., Nets of higher-dimensional cubes it follows that an alternative formulation of the problem is counting spanning trees of the hypercube graph modulo its automorphism group.

@fsl, talking of literature references, A005725 points to an exercise in Comtet's Advanced Combinatorics which gives an integral expression for the coefficient of $x^k$ in $(1 + x + \cdots + x^{r-1})^m$.

Gröbner bases? Or did you specifically mean something specialised to quadratics?

Strengthening a conjecture usually makes it easier to prove, so it is worth noting that empirically the given polynomial for $M = 2m+1$ appears to be $(m+1) (k - 1)^{D(m)}$ (and, FWIW, this holds for $M=1$ too).

@InesInstitoris, I don't know. The closest I've been able to find in about twenty minutes' research is the Hankel continued fractions of this paper, but they should have numerators which are cubic monomials rather than cubic polynomials.

Ah, I got the definition of the Fibonacci-Motzkin paths slightly wrong.

@InesInstitoris, the square consists of one term for the level-0 $\rightarrow$ before moving to level 1 and one term for the level-0 $\rightarrow$ after coming back from level 1.