Q:

Why must basis sets used in computations be non-orthogonal?

While constructing basis sets for calculations using Hartree Fock, the basis set must be non-orthogonal. And because of this, we go through many steps to put the overlap and density matrices in the 'right' basis. (See Modern Quantum Chemistry (Szabo and Ostlund) or David Sherill's notes online). 
My question is: why must this basis be non-orthogonal?

A:

It does not have to be non-orthogonal. The basis set can be arbitrary, and one usually uses the chemically motivated Ansatz with atomic orbitals (LCAO), which just happen to be non-orthogonal. An alternative would be to use a plane wave basis set, which are orthogonal.

A:

The basis doesn't have to be non-orthogonal, but, in practice, it is more efficient to span the desired Hilbert space for each atom with Gaussians than with proper orbitals as Gaussian integrals can be computed relatively quickly. Moreover, if you want to have an atom-centered basis set (which is the way most chemists think about orbitals), these atom-centered basis functions will overlap if they are near each other regardless of what function type you choose, leading to non-orthogonality.

