Orthogonal time frequency space (OTFS) modulation is a recently proposed delay-Doppler (DD) domain communication scheme, which has shown promising performance in general wireless communications, especially over high-mobility channels. In this paper, we investigate DD domain Tomlinson-Harashima precoding (THP) for downlink multiuser multiple-input and multiple-output OTFS (MU-MIMO-OTFS) transmissions. Instead of directly applying THP based on the huge equivalent channel matrix, we propose a simple implementation of THP that does not require any matrix decomposition or inversion. Such a simple implementation is enabled by the DD domain channel property, i.e., different resolvable paths do not share the same delay and Doppler shifts, which makes it possible to pre-cancel all the DD domain interference in a symbol-by-symbol manner. We also study the achievable rate performance for the proposed scheme by leveraging the information-theoretical equivalent models. In particular, we show that the proposed scheme can achieve a near optimal performance in the high signal-to-noise ratio (SNR) regime. More importantly, scaling laws for achievable rates with respect to number of antennas and users are derived, which indicate that the achievable rate increases logarithmically with the number of antennas and linearly with the number of users. Our numerical results align well with our findings and also demonstrate a significant improvement compared to existing MU-MIMO schemes on OTFS and orthogonal frequency-division multiplexing (OFDM).