If $\psi_1(x_1)\psi_1(x_2)$ is antisymmetric (and I understand this is impossible, since the ground state is not degenerate)
The ground state is degenerate, since both particles have the same $n$ (principal) quantum number and thus the same energy. In general, for $N$ particles, the symmetric and antisymmetric wavefunction may be constructed as
\begin{align}\psi_{_S}&\equiv\sqrt{\frac{N_1!\cdots{N}_k!}{N!}}\sum_P\hat{P}\,\phi_{n_1}(\zeta_1)\phi_{n_2}(\zeta_2)\ldots\phi_{n_N}(\zeta_N)\\[0.1in]\psi_{_A}&\equiv\sqrt{\frac{N_1!\cdots{N}_k!}{N!}}\begin{vmatrix}\phi_{n_1}(\zeta_1)&\cdots&\phi_{n_1}(\zeta_N)\\\vdots&&\vdots\\\phi_{n_N}(\zeta_1)&\cdots&\phi_{n_N}(\zeta_N)\end{vmatrix}\end{align}
respectively, where $\zeta_i$ are the internal degrees of freedom and $N_i$ is the degeneracy of the $i$-th set of degenerated particles (for the antisymmetric part, most usually $N_1!\cdots{N}_k!=1$). In your case (given that you can always write the wavefunction as a product of the spatial and spin parts),
$$\psi_{_A}=\begin{vmatrix}\psi_1(x_1)&\psi_1(x_2)\\\psi_1(x_1)&\psi_1(x_2)\end{vmatrix}=0$$
which is why the spatial antisymmetric part is impossible for the ground state. For fermions this is a natural consequence of the Pauli exclusion principle, since you would allow the possibility of two particles being in the same state, given that the spin part would be symmetric.
Now, for the first excited level there is no restriction about considering both the symmetric and antisymmetric parts, in fact you must consider them both. Just as when you must consider the three posibilities from the triplet spin state
$$\chi_{_T}=\begin{cases}\chi_\alpha\\\chi_\beta\\\chi_+\end{cases}$$
you may consider both solutions (4 in total, as you say),
$$\psi=\begin{cases}\frac{1}{\sqrt{2}}\left[\psi_1(x_1)\psi_2(x_2)+\psi_1(x_2)\psi_2(x_1)\right]\chi_-\\\frac{1}{\sqrt{2}}\left[\psi_1(x_1)\psi_2(x_2)-\psi_1(x_2)\psi_2(x_1)\right]\chi_{_T}\end{cases}$$
the thing is that this is a set of possible solutions, just as what you found out for the spin part with the triplet state, the particles may have this or that state, usually when dealing with fermions the only restriction to take care of is Pauli exclusion principle. You may thus consider them all to construct the total wavefunction.
Now, the thing Trimok says,
Note that, mathematically, you can have a total antisymmetric wave function, without having a specific symmetry in the spatial part or in the spin part, for instance : $$\psi_1(x_1)\psi_2(x_2) \alpha(s_1)\beta(s_2) - \psi_2(x_1)\psi_1(x_2) \beta(s_1)\alpha(s_2)$$
could be misleading. This can be seen if you construct the first excited state from the Slater determinant (the general expression for $\psi_{_A}$), say, $n_1=1$, $n_2=2$, $\alpha(1)$, $\beta(2)$, i.e.
$$\mathcal{S}_1\equiv\frac{1}{\sqrt{2}}\begin{vmatrix}\psi_1(x_1)\alpha(1)&\psi_1(x_2)\alpha(2)\\\psi_2(x_1)\beta(1)&\psi_2(x_2)\beta(2)\end{vmatrix}$$
which is the expression given by Trimok, but the thing, again, is that you must consider all possible solutions for this state, meaning that for $n_1=1$, $n_2=2$, you can have
\begin{align}\alpha(1),&\,\alpha(2)\\\beta(1),&\,\beta(2)\\\beta(1),&\,\alpha(2)\\\alpha(1),&\,\beta(2)\end{align}
you can interchange $n_1,\,n_2$ also if you please (there's no new information). For the first, the Slater determinant pops out the antisymmetric spatial solution times $\chi_\alpha$, the second, the antisymmetric spatial part times $\chi_\beta$, but you may take the third,
$$\mathcal{S}_2\equiv\frac{1}{2}\begin{vmatrix}\psi_1(x_1)\beta(1)&\psi_1(x_2)\beta(2)\\\psi_2(x_1)\alpha(1)&\psi_2(x_2)\alpha(2)\end{vmatrix}$$
and the fourth to build the antisymmetric part times $\chi_+$ as $\mathcal{S}_1+\mathcal{S}_2$ and to build the symmetric part times $\chi_-$ as $\mathcal{S}_1-\mathcal{S}_2$. Here both must be taken in count because of the indistinguishability of particles, considering $\mathcal{S}_1$ or $\mathcal{S}_2$ alone is just insufficient. As I showed first, all this is taken care of if you just factor the spatial and spin parts of the wavefunction and treat each one apart, as you were doing.
Yes it does matter whether it's $n_x$ or $n_y$ with the excited state, which means you've got 2 symmetric and 2 anti-symmetric linearly-independent orbital wavefunctions for the first excited state:
$$\psi_{2,1}(r_1)\psi_{1,1}(r_2) \pm \psi_{1,1}(r_1)\psi_{2,1}(r_2)$$
$$\psi_{1,2}(r_1)\psi_{1,1}(r_2) \pm \psi_{1,1}(r_1)\psi_{1,2}(r_2)$$
Taking into account the spins you've got:
$$|0,0\rangle[\psi_{2,1}(r_1)\psi_{1,1}(r_2) + \psi_{1,1}(r_1)\psi_{2,1}(r_2)]$$
$$|1,s_z\rangle[\psi_{2,1}(r_1)\psi_{1,1}(r_2) - \psi_{1,1}(r_1)\psi_{2,1}(r_2)]$$
$$|0,0\rangle[\psi_{1,2}(r_1)\psi_{1,1}(r_2) + \psi_{1,1}(r_1)\psi_{1,2}(r_2)]$$
$$|1,s_z\rangle[\psi_{1,2}(r_1)\psi_{1,1}(r_2) - \psi_{1,1}(r_1)\psi_{1,2}(r_2)]$$
Which is a total of 8 states, all of which are degenerate.
Best Answer
The first step is to realize that the spatial wavefunctions for two non-interacting particles are sums of the form $\psi_{n_1}(x_1)\psi_{n_2}(x_2)$ and so have energy $E=E_{n_1}+E_{n_2}\sim n_1^2+n_2^2$. Thus the lowest energy states have the lowest sums $n_1^2+n_2^2$.
For the spin-1/2 case, the total wavefunction must be fully antisymmetric since spin-1/2 particles are fermions. There are four possible types of spin wavefunctions for two spin-1/2 particles: \begin{align} \chi^1_1&= \chi^{1/2}_+(1)\chi^{1/2}_+(2)\, ,\\ \chi^1_{-1}&= \chi^{1/2}_- (1)\chi^{1/2}_-(2)\, ,\\ \chi^1_0&= \chi^{1/2}_+(1)\chi_-^{1/2}(2)+\chi^{1/2}_-(1)\chi^{1/2}_+(2)\, ,\\ \chi^0_0&=\chi^{1/2}_+(1)\chi_-^{1/2}(2)-\chi^{1/2}_-(1)\chi^{1/2}_+(2) \end{align} The first three are symmetric, while the last is antisymmetric under interchange of the particle label $1\leftrightarrow 2$. Thus, the first three must be combined with a spatially antisymmetric wavefunction, and the last to a spatially symmetric wavefunction, to produce a full antisymmetric total wavefunction.
There are three possible types of spatial wavefunctions. The first two are symmetric wavefunctions of the form $$ \psi^+(x_1,x_2)=\left\{\begin{array}{lol} \psi_{n_1}(x_1)\psi_{n_2}(x_2)+\psi_{n_2}(x_1)\psi_{n_1}(x_2)& \hbox{for } n_1\ne n_2\\ \psi_{n}(x_1)\psi_{n}(x_2)&\hbox{for } n_1=n_2=n \end{array} \right. $$ There is also an antisymmetric combination for $n_1\ne n_2$: $$ \psi^-(x_1,x_2)=\psi_{n_1}(x_1)\psi_{n_2}(x_2)-\psi_{n_1}(x_2)\psi_{n_2}(x_1)\, . $$ which can occur when $n_1\ne n_2$. The permutation symmetry here is on the spatial position $x_1\leftrightarrow x_2$, which results from the interchange particle labels $1\leftrightarrow 2$.
The lowest energy wave function is $\psi_{1}(x_1)\psi_1(x_2)$ and is space-symmetric. Thus, it must be combined with $\chi^0_0$ to produce the fully antisymmetric state with lowest energy.
The next lowest are $\psi_1(x_1)\psi_2(x_2)\pm \psi_1(x_2)\psi_2(x_1)$. The symmetric combination must be multiplied by $\chi^0_0$ to make it fully antisymmetric. The antisymmetric combination can be multiplied by any one of the $\chi^1_m$ to make it fully antisymmetric. This should be enough to produce what you need for the spin-1/2 case.
For the spin-0 case there is only one single particle wavefunction and thus there is only one product spin wavefunction: $\chi^0_0(1)\chi^0_0(2)$; it is fully symmetric under the exchange of particle labels $1\leftrightarrow 2$. It must therefore always be multiplied by a spatially symmetric wavefunction since spin-0 particles are bosons. You can work out the details on this using the three types of spatial wavefunctions given above.