I think an answer can be found in the great paper On the theory of identical particles (1977). I'll write down the argument as I understand it.
Consider 2 identical particles in 3D. The configuration space is $C=(\mathbb{R}^3 \times \mathbb{R}^3)/\mathbb{Z}_2$, where we quotient by $\mathbb{Z}_2$ since the particles are identical. By going to relative coordinates, $C=\mathbb{R}^3 \times (\mathbb{R}^3/\mathbb{Z}_2)$, where $\mathbb{R}^3$ is associated with the center of mass, and where the second is the space of relative coordinates $x\in \mathbb{R}^3$ which are identified with $-x$. This gives a singular point at $x=0$, and we remove it (we keep the particles apart). Thus, the relative part is $(\mathbb{R}^3-{0})/\mathbb{Z}_2=]0,\infty[\times (S^2/\mathbb{Z}_2)$.
First focus on the part $S^2/\mathbb{Z}_2$. As I mentioned above in the question, closed loops in this space fall in two classes. If we think of the sphere $S^2$ before identifying antipodal points, these classes are classes of paths that 1. begin and end on the same point 2. begin and end on antipodal points.
The paths in class 1, viewed in the total space $(\mathbb{R}^3-{0})/\mathbb{Z}_2$, do not encircle the removed point at 0 (or can be continuously deformed to a path that doesn't). The paths in class 2 encircle the origin once (or can be continuously deformed to a path that encircles it once).
The actual argument can now be made using parallel transport of vectors along closed loops in the configuration space. For each $x$ we define a Hilbert space $h_x$ such that $\Psi(x)\in h_x$ is the value of the wave function at $x$. Thus we think of $\Psi(x)$ as a vector that we will parallel transport. Defining a basis $\chi_x$ for $h_x$, we can write $\Psi(x)=\psi(x)\chi_x$. We have a choice of basis, which is a gauge freedom.
We then define a linear, unitary operator $P(x',x) : h_x \to h_{x'}$ which parallel transports vectors of $h_x$ to vectors of $h_{x'}$. It is assumed that it takes the form
$P(x+dx,x)\chi_x = (1+i dx^k b_k (x) ) \chi_{x+dx}$
for infinitesimal parallel displacement. The $b_k$ are functions analogous to the vector potential $A_\mu$, and define gauge-invariant differentiation
$D_k = \partial_k -i b_k (x) $.
In turn, this defines a "curvature" tensor $f_{kl} = i [D_k ,D_l]$ analogous to $F_{\mu \nu}$ in EM or $R_{\mu \nu}$ in GR. The idea, I think, is that parallel transport is trivial along paths which enclose no curvature, but that it can have an effect if the region enclosed is curved.
In our case, $f_{kl}$ is chosen zero for all points in the configuration space, but it is not defined for the excluded point at the origin.
Suppose we consider a path that does not encircle the origin. In that case, the enclosed "curvature" is zero, and we have $P(x,x)\equiv P(x) = 1$. This holds for any path that does not encircle the origin - deforming a path continuously will retain this property. If we consider a path that encircles the origin once, we write $P(x) = \exp[i\zeta]$, since $h_x$ is one-dimensional. Again, if we continuously deform such a path, we cannot change the enclosed "curvature", so the phase will be the same for all such paths.
Then a path corresponding to a "double exchange" yields $P(x) = \exp[2i\zeta]$. Now we continuously deform it back to the trivial path, which we can do without passing through the origin. The "curvature" enclosed for the "double exchange" path is therefore the same as that of the trivial path, namely zero (there should be some nicer way of looking at this). Since for all such paths $P(x)=1$, we have $\exp[2i\zeta]=1$.
If we consider the 2D case, we would get infinitely many classes of paths, since the configuration space would be $\mathbb{R}^2 \times (]0,\infty[\times S^1/\mathbb{Z}_2)$. Yet, the same idea should apply.
