Aligned subdivisions are not that much of a problem actually.
I've tested median cut that can cut at various angles and variants that cut out spheres, and it didn't make much difference.
Choice of the box to cut and the location where you cut is more important. Wu's method is nice in this case because it exactly measures variance of sets after a split, while median cut just estimates.
However, Wu's method needs lookup tables for the whole color (hyper)cube, which requires either large (or for RGBA still ridiculous) amount of RAM or reduced precision, and losing bits of input is way worse than having suboptimal division.
The biggest wins I've found (and implemented in libimagequant) were:
• Voronoi iteration after generating palette (like a final step in LBG quantization). This shifts axis-aligned subdivisions to local optimum (I've tried applying that after every subdivision, but it didn't matter — once in postprocessing is enough!)
• Adjust weights/counts in histogram from quality feedback measured from remapped image. If you're aiming to minimize MSE, then that's not too expensive and it's similar to gradient descent.