We revisit the many-body expansion (MBE) for water-water interactions by examining the effects of the basis set, including those resulting from the basis set superposition error (BSSE) correction, and electron correlation on the various terms for selected sizes of water clusters up to n = 21. The analysis is performed at the second-order Møller-Plesset (MP2) perturbation theory with the family of augmented correlation consistent basis sets up to five zeta quality (aug-cc-pVxZ, x = D, T, Q, 5) for the (H2O)n, n = 7, 10, 13, 16, and 21, clusters for which we report either the complete MBE (for n = 7, 10) or the ones through the 6-body (for n = 13) and the 5-body terms (for n = 16, 21). For the n = 3 and 7 clusters, we also report the analysis at the coupled cluster with single, double, and perturbative triple replacements in order to assess the effects of a higher correlation on the magnitude and percentage of the various MBE terms. Our results suggest that the oscillatory behavior around zero found for the 5-body and larger terms is solely an artifact of the (small) size of the basis set. Indeed, all terms above the 4-body converge monotonically to practically zero upon increasing the size of the basis set toward the complete basis set (CBS) limit. In that respect, the BSSE-corrected 5-body and above terms do not exhibit the oscillatory behavior on either side of zero with the basis set observed for the BSSE-uncorrected terms. In addition, the magnitudes of the 5-body and above terms are accurately reproduced even with the smaller basis set of the series (aug-cc-pVDZ) once the BSSE correction is taken into account. The same level of theory (MP2/aug-cc-pVDZ, BSSE-corrected) also accurately reproduces the MP2/CBS values of the 3- and 4-body terms. The contribution of electron correlation to the 3- and 4-body terms is quite small so that neglecting the correlation contribution in all terms above the 3-body results in an error of the order of 0.1%. The BSSE correction to the largest 2-body term in the MBE was accurately estimated from the function a[1 + erf( - b·R)], which is proportional to the common (overlapping) area between two Gaussian distributions whose centers are separated by R with the constants a and b fitted to the calculated BSSE corrections for the individual 2-body terms of the clusters with each basis set and R is the distance between oxygen atoms. Our results demonstrate that the MBE for water-water interactions converges by the 4-body term since any finite terms above the 4-body are artifacts of the size of the basis set. The MBE can thus be safely truncated at the 4-body term when either a very large basis set is used or BSSE corrections are taken into account even with the smaller aug-cc-pVDZ basis set. We expect these findings to have important consequences in the pursuit of accurate ab initio based many-body molecular dynamics simulations for aqueous systems.