Reverse-time migration of multiples can generate imaging results with wider illumination, higher folds, and broader wavenumber spectra than the conventional migration of primaries. However, the results of the migration of multiples retain heavy crosstalks generated by interactions between unrelated multiples, thereby seriously degrading imaging qualities. To eliminate such crosstalks, we propose a least-squares optimized algorithm of multiples. In this method, different-order water-column multiples and water-bottom-related multiples are extracted using multiple decomposition strategies before migration procedures. The proposed method treats the nth-order water-column multiples as virtual sources for Born modeling to produce the predicted (n+1)th-order water-bottom-related multiples. In each iteration, the gradients are calculated by crosscorrelating the forward-propagated nth-order water-column multiples with the backward-propagated seismic residuals between the observed and predicted (n+1)th-order water-bottom-related multiples. The developed approach is referred to as the least-squares reverse-time migration of water-bottom-related multiples (LSRTM-WM). Numerical experiments on a layered model and the Pluto 1.5 model demonstrate that LSRTM-WM can significantly remove crosstalks and considerably improve spatial resolution.