Molecular orbitals are simply used to construct $N$-electron basis states (Slater determinants).
First, let us have a one-particle Hilbert space $\cal{H}$ of dimension $M$, with $M$ linearly independent vectors
$\{\phi_p\}_{p=1}^M$ forming a basis on it. In non-relativistic one-electron QM, $\cal{H}=L^2(\mathbb{R}^3)\otimes\mathbb{C}^2$ (scattering states now not considered for simplicity), so in principle, $M\rightarrow\infty$ should be taken; we nevertheless keep $M$ finite to properly see the dimensions (also, $M$ is obviously finite in all numerical calculations). For the sake of simplicity, let us further take the basis to be orthonormal:
$$\langle\phi_p|\phi_q\rangle=\delta_{pq} \ , $$
although this is not necessary. Since this is a basis, any one-particle function can be expanded in it as
$$
|\psi\rangle=\sum_{p=1}^Mc_p|\phi_p\rangle \ .
$$
The important thing to realize is that constructing the antisymmetrized products of $N$ such vectors in all ${M}\choose{N}$ possible ways will leave you with an orthonormal basis on the antisymmetric subspace of the $N$-particle Hilbert space. For example, for $N=2$, the antisymmetrized products
$$
|\Phi_{pq}\rangle=\frac{1}{\sqrt{2}}\left[
|\phi_p\rangle\otimes|\phi_q\rangle-|\phi_q\rangle\otimes|\phi_p\rangle
\right]
$$
for $p,q=1,...,M$, $p<q$, form a basis in the antisymmetric subspace of the two-particle Hilbert space $\cal{H}\otimes\cal{H}$ (which is the relevant subspace for two electrons). Any two-electron function can be expanded in this ${{M}\choose{2}}=M(M-1)/2$ -dimensional basis:
$$
|\Psi\rangle=\sum_{\substack{p,q=1 \\ p<q}}^Mc_{pq}|\Phi_{pq}\rangle \ ,
$$
orthonormality being understood as
$$
\langle\Phi_{rs}|\Phi_{pq}\rangle=\delta_{pr}\delta_{qs}-\delta_{ps}\delta_{qr} \ .
$$
The idea is the same for the general $N$-electron case, the antisymmetrized products (Slater determinants)
$$
|\Phi_{p_1p_2...p_N}\rangle=
\frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}}
\left[|\phi_{p_1}\rangle\otimes|\phi_{p_2}\rangle\otimes...\otimes|\phi_{p_N}\rangle\right]
$$
can be used in
$$
|\Psi\rangle=\sum_{\substack{p_1,p_2,...,p_N=1 \\ p_1<p_2<...<p_N}}^Mc_{p_1p_2...p_N}|\Phi_{p_1p_2...p_N}\rangle \ .
$$
Such expansions are used extensively in the many-body theory.
In the most basic version of Full-CI, these ${M}\choose{N}$ Slater determinants are used to construct the matrix representation of the Hamiltonian, and the energies and coefficients $c_{p_1p_2...p_N}$ of the eigenstates are obtained from the diagonalization of this matrix.
Note that, apart from orthonormality, I assumed nothing about the one-particle basis. These orbitals are usually found from e.g. Hartree-Fock or Kohn-Sham calculations. On the level of theory, LCAO does not have much to do with any of this, it is just often (almost always) useful to expand molecular orbitals over atom-centered functions:
$$
|\phi_p\rangle=\sum_{\mu=1}^Md_{\mu p}|\chi_\mu\rangle \ .
$$
This is practical from a computational/interpretational point of view, as it is often easier to see e.g. bonding/antibonding orbitals in this way.
But the expansion over atomic functions is by no means necessary, you can do quantum chemistry in e.g. a plane wave basis, too.
Update
I used bra-ket notations above in order to properly show the tensor product structure. The (probably) better known form of these functions is recovered by projection with the formal coordinate eigenstates $|\vec{r}\rangle$. The one-electron wave functions (orbitals) then read
$$
\phi_{p}(\vec{r})=
\langle\vec{r}|\phi_p\rangle=
\left(
\begin{matrix}
\phi_{p\uparrow}(\vec{r}) \\
\phi_{p\downarrow}(\vec{r})
\end{matrix}
\right) \ ,
$$
having two components due to their spinor character (this is the $\mathbb{C}^2$ part of ${\cal{H}}$). The components $\phi_{p\uparrow}(\vec{r})$ and $\phi_{p\downarrow}(\vec{r})$ are different in the general case; one usually wants orbitals to be eigenfunctions of $\hat{s}_z$, which forces one of the components to be zero.
The $N$-electron Slater determinants are obtained similarly after a projection with
$$
|\vec{r}_1,\vec{r}_2,...,\vec{r}_N\rangle=|\vec{r}_1\rangle\otimes|\vec{r}_2\rangle\otimes...\otimes|\vec{r}_N\rangle \ ,
$$
leading to
$$
\Phi_{p_1p_2...p_N}(\vec{r}_1,\vec{r}_2,...,\vec{r}_N)=
\frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}}
\left[\phi_{p_1}(\vec{r}_1)\otimes\phi_{p_2}(\vec{r}_2)\otimes...\otimes\phi_{p_N}(\vec{r}_N)\right] \ ,
$$
which shows that non-relativistic $N$-electron states have $2^N$ spinor components.
Note that ${\cal{P}}$ permutes the orbital indices, not the coordinate indices.
When doing formal manipulations with wave functions, it is often useful to introduce a combined spatial+spin coordinate eigenstate $|x\rangle=|\vec{r},\sigma\rangle$ for $\sigma=\uparrow,\downarrow$, which implements a further projection onto the upper or lower spinor component:
$$
\langle\vec{r},\uparrow\hspace{-0.1cm}|\phi_p\rangle=\phi_{p\uparrow}(\vec{r})
\ \ \ , \ \ \
\langle\vec{r},\downarrow\hspace{-0.1cm}|\phi_p\rangle=\phi_{p\downarrow}(\vec{r}) \ .
$$
Projecting $|\Phi_{p_1p_2...p_N}\rangle$ with
$$
|x_1,x_2,...,x_N\rangle=|x_1\rangle\otimes|x_2\rangle\otimes...\otimes|x_N\rangle
$$
will similarly pick out a single one from the $2^N$ spinor components:
$$
\Phi_{p_1p_2...p_N}(x_1,x_2,...,x_N)=
\frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}}
\left[\phi_{p_1}(x_1)\phi_{p_2}(x_2)...\phi_{p_N}(x_N)\right] \ .
$$
In this way, it is often easier to talk about the antisymmetry of the electronic wave function under simultaneous interchanges in coordinate and spin space:
$$
\Psi(x_1,x_2)=-\Psi(x_2,x_1) \ .
$$
