We derive a new multistage sphere decoding (MSD) algorithm, a generalization of the conventional sphere decoder (SD). This new MSD exploits that many higher-order signal constellations can naturally be decomposed into several lower-order constellations. We develop a 2-stage SD for a 16-ary quadrature amplitude modulation (16QAM) MIMO system by decomposing 16QAM into two 4QAM constellations. The first stage generates a list of 4QAM vectors. For each of these, the second stage computes an optimal 4QAM vector. In the low SNR region, our MSD performs close to the original (single-stage) SD, but has lower complexity. In the high SNR region, our MSD is not suitable for reaching near ML performance.