Why are complex basis functions not used in electronic structure calculations?
The basic answer is rather simple: In the non-relativistic limit it is just not necessary. There is no complex operator and hence nothing which could change the complex component of your basis functions. Using complex basis functions would only double the required work (i.e., computer time). There are very efficient schemes in order to calculate the standard integrals (McMurchie-Davidson, Obara-Saika, ...). I don't think that the representability of MOs is a reason for using real basis functions, since, strictly speaking, you can interprete only the square of the MOs and not their amplitude.

However, there are situations, where your MOs can be complex. This can happen when you have a complex operator in your equation. An example are the Pauli spin matrices of which $\sigma_y$ has only an imaginary compound. They occur, e.g., the Dirac equation (see here one formulation where they explicitly state the incorporation of the Pauli matrices). Complex MOs can be obtained either by using complex basis functions and real MO coefficients or by using real basis functions with complex coefficients. Starting from a non-relativistic implementation it might be simpler just to keep your basis functions real and to choose the MO coefficients complex.