How do we measure distances in the FLRW metric? We are in the flat FLRW metric in Cartesian comoving coordinates. The metric is expressed as:
$$ds^2 = d\tau^2 + a(\tau)^2\big(dx^2 + dy^2 + dz^2\big)$$
The fact that the "universe is expanding" is justified by saying that if I have a galaxy at $(\tau, 0, 0, 0)$ and one at $(\tau, L, 0, 0)$ and I measure the distance $R(\tau)$ along the curve $\sigma_\tau : x \mapsto (\tau, x, 0, 0)$ from $x=0$ to $x=L$, I find that this distance $R(\tau)$ is indeed a function of $\tau$.
However, isn't the distance between two spacelike-separated events defined as the lengths of the geodesic connecting them? 
This curve is not a geodesic. I was definitely surprised. 
To see this, let's look ay the geodesic equations
$$ \frac{d^2\gamma^{\nu}}{d\lambda^2} = -\Gamma^{~\nu}_{\mu\rho}\frac{d\gamma^{\mu}}{d\lambda}\frac{d\gamma^{\rho}}{d\lambda}$$
with initial conditions 
\begin{aligned}
\gamma^\mu(0) &= (\tau, 0, 0, 0)\\\\
\left.\frac{d\gamma^{\mu}}{d\lambda}\right|_{\lambda=0} &= (0, 1, 0, 0)
\end{aligned}
Remember that the only non-zero Christoffel symbols are 
$\Gamma^{~\tau}_{ii} = a\dot a$ and
$\Gamma^{~i}_{i\tau} = \Gamma^{~i}_{\tau i} = \frac{\dot a}{a}$ for $i = x, y, z$, and where the dot represents differentiation by $\tau$. The geodesic equations become
\begin{aligned}
\frac{d^2\gamma^{\tau}}{d\lambda^2} &= -a\dot a \left(\frac{d\gamma^{x}}{d\lambda}\right)^2 \\ \\
\frac{d^2\gamma^{x}}{d\lambda^2} &= -\frac{\dot a}{a}\frac{d\gamma^{x}}{d\lambda}\frac{d\gamma^{\tau}}{d\lambda}\\
\end{aligned}
In particular, we see that for $\lambda = 0$, 
$$\frac{d^2\gamma^{\tau}}{d\lambda^2} = -a\dot a$$ from which we learn two things:


*

*our first curve $\sigma_\tau$ is not a geodesic

*spacelike geodesics in the FLRW metric do not tay on surfaces of constant cosmological time $\tau$


Which brings me to my question:


*

*What is the relation between the distance $R(\tau)$ and the proper distance (the one measured along a spacelike geodesic?)


Or, perhaps in a different formulation:


*

*How does one make sense physically of $R(\tau)$? and of the proper distance?

 A: From my calculation you are correct that it is not a geodesic [Update: see reference below]. Consider a general (not necessarily spatially flat) Friedmann-Lemaitre-Robertson-Walker spacetime. I use the "hyperspherical" coordinates, with line element
$$ds^2=-dt^2+a(t)^2\big(d\chi^2+S_k^2(d\theta^2+\sin^2\theta\,d\phi^2)\big)$$
Consider the curve $x^\mu=(t_0,\chi,\theta_0,\phi_0)$, so only the radial coordinate $\chi$ varies and we also take $\chi$ as the parameter. The tangent vector in our coordinates is simply $(0,1,0,0)$. Now the geodesic equation written in the original question is only valid for an affine parameter, but in my experience many sources fail to qualify this clearly. Eric Poisson gives an example in FLRW with non-affine parameter (27:30 to 56:00; you may wish to check this to see how relevant it is, and summarise for the rest of us). So let's follow Poisson's textbook A Relativist's Toolkit, $\S1.3$.
Write $\mathbf u$ for our spatial tangent vector above, then as I understand $\nabla_{\mathbf u}{\mathbf u}\propto\mathbf u$ for a geodesic, and in particular the LHS is zero if the parameter is affine. However I get the vector
$$(a\dot a,0,0,0)$$
for the LHS, which is certainly not parallel to $\mathbf u$. This concurs with the conclusion of the original poster. I did not find mention of this counter-intuitive result after a quick look in MTW, Griffiths & Podolsky, and Plebanski & Krasinski. However, we would certainly still expect that within the 3-dimensional hypersurface the curve would be a geodesic!
Now speaking of 3-dimensional hypersurfaces, in an arbitrary spacetime we can define these as being orthogonal to the worldlines of certain preferred observers, in cosmology take the motion of matter as averaged over some local region but see e.g. Ellis+ Relativistic Cosmology, $\S4.6.2$ for potential issues. But in FLRW it is straightforward, take the observers who measure an isotropic CMB. Once we have such a 3+1-splitting of spacetime, measure spatial distance along these hypersurfaces. From a "physics" point of view perhaps this is more practical than sending out spatial geodesics from an individual observer's worldline?
Update (August 2018): Ellis writes in the 2012 paper "The evolving block universe and the meshing together of times":

...FLRW spacetimes (where the surfaces of homogeneity are generically not simultaneous according to the radar definition).

Simultaneity via radar is the same as Einstein simultaneity, both yield the spacelike hypersurface orthogonal to the given observer's worldline; at least if done locally. So if Ellis means locally, then he concurs. See also the reference cited therein, which is Ellis & Matravers 1985, "Spatial homogeneity and the size of the universe."
A: I'll answer the second bulleted version of your question, and argue that answering the first version isn't worth the effort.
Here's a simple heuristic illustration of what's going on. In a spatially closed, Euclidean-time universe that's expanding at a constant rate, space"time" is a cone. The cosmological time (the time coordinate with respect to which this geometry is homogeneous at any fixed time) is the straight-line distance from the apex. If you unroll the cone and lay it flat, the circles of constant time are circular arcs and the distance between points on an arc is obviously shorter in the background spacetime than along the arc.
Here's a numerically precise example: in the $\Omega=0$ limit of FLRW cosmology, in certain units, $a(\tau) = \tau$ and the metric is $ds^2 = d\tau^2 - \tau^2 (d\chi^2 + {\sinh^2}\chi \, (d\theta^2 + {\sin^2}\theta\,d\phi^2))$. This is just a part of Minkowski space (specifically the future light cone of some point) as you can see by doing the substitution $t = \tau\cosh\chi$, $r = \tau\sinh\chi$, $\theta=\theta$, $\phi=\phi$. The surfaces of constant $\tau$ are spheres (of timelike radius) centered at the origin.
It's true generally that FLRW coordinates are analogous to polar coordinates, surfaces of constant cosmological time are analogous to spheres, and spacelike geodesics between points on those surfaces are analogous to chords of spheres.
I would expect the relationship between the constant-$\tau$ distance and the "chord distance" to be complicated in general since $a(\tau)$ is. At least in the $\Omega{=}0$ case it could be worked out easily, but I don't think it would provide any physical insight, because spacelike distances in general have no physical significance. Even, e.g., the proton radius is an abstraction of the consequences of light-speed-limited processes that only see spacetime at much smaller scales; it's useful to human physicists but doesn't enter into the physical laws in any way. Cosmology has a large-scale "starburst" shape simply because the center of it did for some reason and the rules governing the rest of it are rotationally symmetric. The presence of large concentric circles in the pattern is not coincidence, but it doesn't really mean anything.
