The question is confusing two things: the linear combination of atomic orbitals (LCAO) scheme where the molecular orbitals are expanded as a sum of the atomic orbitals
\begin{equation}
\psi_{i\sigma} ({\bf r}) = \sum_{\alpha} C_{\alpha i} \chi_\alpha ({\bf r}) \label{eq:lcao}
\end{equation}
and the actual definition for the atomic basis functions, which in the case of contracted Gaussian basis functions is
\begin{equation}
\chi_\alpha ({\bf r}) = r^{l_\alpha} \left[\sum_{k} d_{k \alpha} \exp(-\zeta_kr^2)\right] Y_{l_\alpha}^{m_\alpha}(\hat{\bf r})
\end{equation}
where I have assumed a pure spherical basis for ease of notation. The LCAO scheme is not limited to Gaussian basis sets; the linear expansion ansatz is used in every basis set method, regardless of the nature of the basis functions (e.g. Gaussian or Slater-type orbital, numerical atomic orbital, plane waves, finite elements, etc). I have recently written two open access review articles on these topics, see Int. J. Quantum Chem. 119, e25968 (2019) for discussion on various kinds of basis functions and Molecules 25, 1218 (2020) for the self-consistent field algorithm.
It is true that the contracted functions in many basis sets such as the cc-pVXZ basis sets represent occupied atomic orbitals in Hartree-Fock calculations. However, this is not true for all basis sets. As programs tend to aim for generality, they typically aim to implement initial guesses that are not basis set dependent, such as the superposition of atomic densities - which is equivalent to occupying the contracted functions in certain basis sets but also works for basis sets where the individual basis functions do not look like occupied orbitals. Refer to J. Chem. Theory Comput. 15, 1593 (2019) for further discussion on initial guesses.
I would like to note that unlike stated in the comments, the matrix ${\bf C}$ is not square in general; this is unfortunately still a very common misconception in computational chemistry. There are often linear dependencies in the basis set (e.g. a basis function on one atom can be also well-described by the sum of basis functions on the neighboring atoms), and removing them removes some columns from ${\bf C}$. Stereotypical cases of this happening are when you use diffuse functions that are needed to describe long-range interactions and weakly bound electrons, or when two nuclei are close together. However, even such pathological linear dependencies can be cured and Coulomb barriers to nuclear fusion modeled when suitable numerical algorithms are chosen to pick a suitable sub-basis (see the links). The choice of the linearly independent basis is also discussed in the overview on SCF theory linked above.
