This paper investigates multi-stream downlink precoding for massive multiple-input multiple-output low-Earthorbit satellite (SAT) communication systems. We adopt a delay and Doppler precompensation approach to achieve coherent transmission. Under this setting, we formulate a signal transmission model that incorporates the near-independent properties of inter-SAT interference and compensation errors. We then demonstrate that moving beyond single-stream transmission requires both multi-SAT cooperation and multi-antenna UTs. Based on this configuration and the established signal transmission model, we derive the first- and second-order statistical channel characteristics and utilize them to design locally optimal precoding algorithms for both total power constraint (TPC) and per-antenna power constraint (PAPC) conditions, which rely only on statistical channel state information (sCSI). In particular, the designed PAPC algorithm achieves linear complexity with respect to the number of antennas on the cooperative SATs. To reduce the computational complexity of the locally optimal precoder under TPC, we propose a low-complexity and robust precoding scheme optimized for both minimum mean squared error and sum-rate maximization objectives. Using majorization theory, we also provide a rigorous theoretical analysis of the optimal precoding structure under TPC. Moreover, the Lanczos algorithm is adopted to further reduce the complexity of the proposed robust designs. Simulation results show that when each SAT is equipped with a sufficiently large number of antennas, the proposed sCSI-based designs achieve performance comparable to that of instantaneous CSI-based designs.