The case for Gaussians is purely computational. The product of two Gaussians on different centers is a single Gaussian on a third (the product theorem), so every multi-center one- and two-electron integral collapses to a standard closed form. Slater functions have the correct shape — finite cusp, exponential tail — but obey no such identity, so their multi-center integrals must be done numerically and are far slower. In a calculation dominated by the $O(M^4)$ two-electron integrals, that speed gap decides everything.

The shape is then bought back in pieces. A contracted basis function is a fixed sum of primitives — STO-3G uses three Gaussians per orbital, fitted to mimic one Slater function. More primitives, or more contracted functions (double-zeta, triple-zeta), approach the true shape as closely as you are willing to pay for. A finite Gaussian basis never makes the cusp exactly, but it gets close enough that the residual is smaller than the method error you are already accepting. Wrong shape, right answer, fast integrals — the working compromise of the whole field.