This answer has two parts. The first part gives partial answers, and the second part tries to complete the picture using a different formulation.
Part 1
(1.) ... Can the state of a single particle be symmetric or anti-symmetric?
No. As the OP suspected, the "fermion" or "boson" character of particles is relevant only when multiple particles of the same species are present. Note that the context of the question is nonrelativistic quantum mechanics. In relativistic QFT, the "number of particles" is typically not well-defined.
(2.) What does spin even mean in the context of symmetrized Fock-spaces?
See Part 2.
(3.) What is the right single particle Hilbert space when I want to consider e.g. a gas of $N$ electrons?
For $N$ non-relativistic fermions, see Part 2.
(4.) In opposite to (3.), what would be the right single particle space if I want to describe N bosons (what are bosons?)? E.g. photons.
For $N$ non-relativistic bosons (which implies massive bosons), see Part 2.
For photons, the story is different. Photons are never non-relativistic. The tensor product construction is appropriate for non-relativistic models, but not for relativistic models. For the latter, quantum field theory is more appropriate. Quantum field theory is not formulated in terms of particles at all; it is formulated in terms of fields. Particles, like photons, are phenomena that the theory predicts.
The construction shown in Part 2 may be regarded as non-relativistic quantum field theory, where the connection between fields and particles is more straightforward.
(5.) If spin is just an eigenvalue of an operator on Fock-space, how do I know that the differences of eigenvalues in bosonic systems are integer and in fermionic systems are half-integer?
This connection does not come from non-relativistic QM, where it is merely enforced by hand. The connection comes from relativistic QFT, where it is called the spin-statistics theorem. For more information and some references, see this post:
https://physics.stackexchange.com/a/441790.
(6.) What is the right Hilbert space for a single boson?
For a single non-relativistic boson with spin $0$, $L^2(\mathbb{R}^3)$ is sufficient. For higher-spin bosons, see Part 2.
Part 2
Up to isomorphism, there is only one infinite-dimensional separable Hilbert space. When we construct that Hilbert space in a particular way, what we're really doing is setting it up so that local observables can be expressed in a relatively simple way. The construction described here takes that idea a step further, by first constructing an abstract algebra that contains all the local observables, and then using this algebra to construct a Hilbert-space representation. After this is done and the Hamitonian has been specified, the connection with spin will (hopefully) be made clear.
Some caveats:
This construction is specific to nonrelativistic models (so it excludes photons), but a similar observables-first approach is also useful in relativistic QFT.
The approach can be used for an arbitrary number of different species
of fermions and bosons, each of arbitrary spin. Fermions will be used here as an example, with a brief comment about how to modify things for bosons.
The algebra will be specified in terms of
generators and relations. This will (initially) be an abstract
*-algebra, but I'll still call its elements "operators".
Initially, instead of "species" and "spin components", the
formulation uses just one index to indicate different fermion
operators. The distinction between "species" and "spin components" is introduced through the Hamiltonian that governs time-evolution.
To avoid technical complications,
continuous space will be replaced by a
very-fine-but-discrete spatial lattice $L$ of
very-large-but-finite extent, so that the total
number of lattice sites is finite.
The "continuum limit" will be treated heuristically
because we don't actually need continuous-space concepts for anything
other than motivation.
Some continuum-inspired notation will be used, like
\begin{equation}
\delta^D(x-y) \equiv \epsilon^{-D}\delta_{xy}
\hskip2cm
\int d^Dx\ f(x)\equiv \epsilon^D\sum_x f(x),
\tag{1}
\end{equation}
where $\epsilon$ is the lattice spacing, $D$ is the number
of spatial dimensions ($D=3$), and $x,y\in L$ are lattice sites.
Here's the construction of the algebra:
With each lattice site $x$,
associate operators $a_n(x)$ with adjoints $a_n^*(x)$
satisfying
\begin{align*}
a_n(x)a_k^*(y)+a_k^*(y)a_n(x) &=\delta_{nk}\delta^D(x-y)
\tag{2}
\\
a_n(x) a_k(y)+a_k(y)a_n(x) &= 0.
\tag{3}
\end{align*}
Together with the usual properties like associativity
and linearity, the relations (2)-(3) define a *-algebra
generated by the $a$s. Equations (2)-(3) are for fermions.
For bosons, replace "+" with "$-$" on the left-hand side of both equations.
If fermion and boson operators are both included,
they commute with each other.
To construct a Hilbert-space representation
of this algebra, let $|0\rangle\in{\cal H}$
be a vector (which will end up being
the vacuum vector) that satisfies
\begin{equation}
a_n(x)|0\rangle = 0
\hskip2cm
\langle 0|0\rangle = 1
\tag{4}
\end{equation}
for all $n,x$.
All other vectors in the representation are linear combinations of
vectors obtained from $|0\rangle$ by applying sums of products
of the adjoint operators $a_n^*(x)$. The inner products
between these vectors may be inferred from equations (2)-(4).
Presumably this inner product space can be completed to define
a Hilbert space on which the operators $a_n(x)$ are densely defined
(they are unbounded), but I won't try to prove this.
To hide the fact that we're working on a discrete spatial lattice
instead of in continuous space,
we can define the smeared field operators
\begin{equation}
a_n(f)\equiv \int d^Dx\ f(x)a_n(x)
\hskip2cm
a_n^*(f)\equiv (a_n(f))^*.
\tag{5}
\end{equation}
(Recall the definition (1) of the "integral.") These satisfy
\begin{align*}
a_n(f) a_k^*(g)+a_k^*(g)a_n(f) &= \delta_{nk}\int d^Dx\ g^*(x)f(x)
\tag{6}
\\
a_n(f) a_k(g)+a_k(g)a_n(f) &= 0.
\tag{7}
\end{align*}
If we take the smearing functions $f,g$ to be "smooth" enough
(slowly-varying compared to the lattice spacing $\epsilon$),
then we might as well be working in continuous space.
So far, we have a *-algebra of operators represented on an inner product space,
but nothing has been said about the physical significance of these
operators. In particular, nothing has been said about spin, which is what the question is all about.
Also, nothing has been said yet about time. Let's fix this first.
An example of a Hamiltonian operator $H$ will be given below.
Given such an operator,
we can construct a one-parameter family of
unitary operators $U(t)=\exp(-iHt)$,
and we can define time-dependent versions of the field operators like this:
\begin{equation}
a_n(f,t)\equiv U(-t)a_n(f)U(t).
\tag{8}
\end{equation}
Since this is a unitary transformation, equation (6) implies
\begin{align*}
a_n(f,t) a_k^*(g,t)+a_k^*(g,t)a_n(f,t) &= \delta_{nk}\int d^Dx\ g^*(x)f(x)
\tag{9}
\end{align*}
and similarly for equation (7).
These are the equal-time (anti)commutation relations.
The operator $H$ will be chosen
so that $H|0\rangle=0$, and then equation (4) implies
\begin{equation}
a_n(f,t)|0\rangle = 0.
\tag{10}
\end{equation}
Now, here is the beginning of a physical interpretation:
Given a spatial region $R\subset L$, the subalgebra
generated by the operators $a_n(f,t)$ and their adjoints
represent operators that are "localized" in $R$ at time $t$
if the supports of the smearing functions $f$ are restricted to $R$.
Within this subalgebra, operators like $a_n^*(f)a_n(g)$
that are invariant under
\begin{align*}
a_n(f)\rightarrow e^{i\theta}a_n(f)\hskip2cm\text{for all }n,f
\tag{11}
\end{align*}
for all $\theta\in \mathbb{R}$ qualify as observables
localized in $R$ at time $t$.
These are local observables. Note that observables
localized in non-overlapping spatial regions commute with each other,
because they necessarily involve products of an even number of $a$s and $a^*$s.
Applying the operator $a_n^*(f,t)$ to any state-vector
creates a particle with wavefunction $f$ at time $t$,
and then equations (2)-(4) imply
that applying the operator $a_n(f,t)$ to any state-vector
annihilates a particle, or gives the zero-vector
if there was no compatible particle present to begin with.
In a single-particle state, or in a state involving only
particles with different indices, the fermion/boson character
is irrelevant; but in a generic multiparticle state, it matters.
All that remains is to specify the Hamiltonian $H$, and this is finally
where spin enters the picture.
A typical Hamiltonian (again for a non-relativistic model) is
\begin{align*}
H &= \int d^Dx\
\left(
\sum_{j,n}\frac{
(\nabla_j a_n)^*(x)\,(\nabla_j a_n)(x)
}{2m_n}
+
\sum_{j,n,k} a_n^*(x) \sigma_{nk}^j B^j(x)a_k(x)\right)
\\
&+\text{terms involving 4 or more $a/a^*$s}
\tag{12}
\end{align*}
where $\nabla_j$ is a finite-difference version of the gradient along
the $j$-th spatial direction.
The unwritten terms
are all chosen so that $H$ is a self-adjoint observable
(recall the requirement (11)) with all annihilation operators $a$
written to the right of all creation operators $a^*$ so that $H|0\rangle=0$.
They represent mutual interactions between the particles.
For now, the quantities $\sigma^j_{nk}$ are just unspecified numeric coefficients.
The functions $B^j(x)$ represent the influence of an external (prescribed)
magnetic field, which will be used to motivate the spin-interpretation.
Even though we defined everything using a discrete lattice
instead of continuous space, space might as well be continuous
if we consider only state-vectors obtained from $|0\rangle$
by applying $a_n^*(f)$s with sufficiently "smooth" smearing functions $f$,
as defined above.
With that in mind, we can pretend (for motivation only)
that the integral in (12) really is an integral over $x$
rather than being secretly a sum over $x$.
Now, finally, here is how spin enters the picture.
If we choose the coefficients $\sigma^j_{nk}$ so that the Hamiltonian
is invariant under a rotation of $(B^1,B^2,...)$ combined
with a transformation $a_n\mapsto \sum_k U_{nk}a_k$ that
leaves the Hamiltonian invariant (in the continuum limit),
then we can use this rotational symmetry to justify interpreting
the different $a_n$s that are mixed with each other by this
transformation as different spin-components of a single species
of fermion.
Different groups of $a_n$s that are not mixed with each other
by this transformation are regarded as different species,
such as electrons versus protons.
Within each group, the mass-parameters $m_n$ be the
same in order for the first term in the Hamiltonian (the kinetic term)
to be invariant.
In the simplest case where the $a_n$s are mixed with each other
only within pairs, with no mixing between the pairs,
we have a system of spin-1/2 fermions, and the coefficients
$\sigma^j_{nk}$ are a representation of the Pauli matrices.
This is one way to define the Pauli matrices, or at least to motivate the their definition.
To relate this formulation to the formulation used in the question,
the key is to notice that the number of particles
(the number of $a_n^*$s applied to the vacuum state)
is invariant under time-evolution because every term in the Hamiltonian has an equal number of $a$s and $a^*$s. As a result, the whole model decomposes
into separate sub-models, one for each different number of particles,
whose state-vectors do not mix with each other.
(This is a special property of strictly non-relativistic models.
In relativistic models, the number of particles is not conserved.)
Each one of these sub-models can be reformulated by expressing
everything — including the effect of the Hamiltonian —
in terms of the smearing functions $f$, which then play the
role of single-particle wavefunctions, multiplied together and
automatically appropriately (anti)symmetrized thanks to
equations (2)-(3). I won't work that out in detail here, because this
post is already obnoxiously long.
