The revolutionary technology of \emph{Stacked Intelligent Metasurfaces (SIM)} has been recently shown to be capable of carrying out advanced signal processing directly in the native electromagnetic (EM) wave domain. An SIM is fabricated by a sophisticated amalgam of multiple stacked metasurface layers, which may outperform its single-layer metasurface counterparts, such as reconfigurable intelligent surfaces (RISd) and metasurface lenses. We harness this new SIM concept for implementing efficient holographic multiple-input multiple-output (HMIMO) communications that dot require excessive radio-frequency (RF) chains, which constitutes a substantial benefit compared to existing implementations. We first present an HMIMO communication system based on a pair of SIMs at the transmitter (TX) and receiver (RX), respectively. In sharp contrast to the conventional MIMO designs, the considered SIMs are capable of automatically accomplishing transmit precoding and receiver combining, as the EM waves propagate through them. As such, each information data stream can be directly radiated and recovered from the corresponding transmit and receive ports. Secondly, we formulate the problem of minimizing the error between the actual end-to-end SIMs'parametrized channel matrix and the target diagonal one, with the latter representing a flawless interference-free system of parallel subchannels. This is achieved by jointly optimizing the phase shifts associated with all the metasurface layers of both the TX-SIM and RX-SIM. We then design a gradient descent algorithm to solve the resultant non-convex problem. Furthermore, we theoretically analyze the HMIMO channel capacity bound and provide some useful fundamental insights. Extensive simulation results are provided for characterizing our SIM-based HMIMO system, quantifying its substantial performance benefits.