So it turns out my earlier intuition was incorrect:

Theorem. Let $X, Y$ be real Hilbert spaces, with the dimension of $X$ at least two, and let $f: X \to Y$ be a function such that $\|f(x)-f(y)\|=\|x-y\|$ whenever $\|x-y\|$ is a natural number. (We do not assume $f$ to be continuous.) Then $f$ is an isometry.

Proof By passing to an affine subplane of $X$, we may assume without loss of generality that $X$ is two-dimensional, so we will write $X = {\bf R}^2$.

If $e$ is a unit vector in $X$, $x$ is an element of $X$ and $n,m$ are natural numbers then $f(x), f(x+ne), f(x+(n+m)e)$ form a degenerate triangle with side lengths $n,m,n+m$, which implies that $f(x+ne) = f(x) + nu$ and $f(x+(n+m)e) = f(x) + (n+m) u$ for some unit vector $u$. In particular, $f$ maps any isometric copy of ${\bf Z}$ to another isometric copy of ${\bf Z}$.

Next, if $e, f$ are orthogonal unit vectors in $X$, then $f(x), f(x+3e), f(x+4f)$ form a triangle with side lengths $3,4,5$ and in particular $f(x+3e)-f(x)$ is orthogonal to $f(x+4f)-f(x)$. From this and the previous claim we see that $f$ maps any isometric copy of ${\bf Z}^2$ to another isometric copy of ${\bf Z}^2$.

By the preceding discussion we have $$f(x,y)-f(0,0) =x (f(1,0)-f(0,0)) + y (f(0,1)-f(0,0)) \qquad (1)$$ whenever $x,y$ are integers, but also $f(n e)-f(0,0) = n (f(e)-f(0,0))$ for any unit vector $n$. In particular, for any Pythagorean triple $a^2+b^2=c^2$, we have (1) when $(x,y)$ is an integer multiple of $(\frac{a}{c},\frac{b}{c})$; conjugating by translation we then see that (1) also holds when $(x,y)$ is an integer multiple of $(\frac{a}{c},\frac{b}{c})$ plus an element of ${\bf Z}^2$. Iterating this we see that (1) holds whenever $(x,y)$ is an integer linear combination of Pythagorean fractions $(\frac{a}{c}, \frac{b}{c})$. In particular, (1) holds on a dense subset $D$ of ${\bf R}^2$.

This is already enough to get isometry in the continuous case. Now we remove the continuity hypothesis. Observe that if $\|x-y\| \leq 2$, then $x$ can be reached from $y$ by a pair of moves of unit length, and hence by the triangle inequality $\|f(x)-f(y) \| \leq 2$.

Now let us normalise $f,X,Y$ so that ${\bf R}^2 \subset Y$, $f(0,0)= (0,0)$, $f(1,0)=(1,0)$, and $f(0,1)=(0,1)$. Then by (1) we see that $f$ is the identity on $D$. For any $p \in {\bf R}^2$ and $\varepsilon>0$, we can find elements $q, r$ of $D$ within $\varepsilon$ of $p + (2,0)$ and $p-(2,0)$ that lie within $2$ of $p$, and thus by the preceding we see that $f(p)$ lies within $2+\varepsilon$ of $p+(2,0)$ and $p-(2,0)$. Sending $\varepsilon$ to zero we see that $f(p)=0$, and the claim follows.

So it turns out my earlier intuition was incorrect, and one can leverage the order properties of ${\bf R}$ to show:

Theorem. Let $X, Y$ be real Hilbert spaces, with the dimension of $X$ at least two, and let $f: X \to Y$ be a function such that $\|f(x)-f(y)\|=\|x-y\|$ whenever $\|x-y\|$ is a natural number. (We do not assume $f$ to be continuous.) Then $f$ is an isometry.

Proof By passing to an affine subplane of $X$, we may assume without loss of generality that $X$ is two-dimensional, so we will write $X = {\bf R}^2$.

If $e$ is a unit vector in $X$, $x$ is an element of $X$ and $n,m$ are natural numbers then $f(x), f(x+ne), f(x+(n+m)e)$ form a degenerate triangle with side lengths $n,m,n+m$, which implies that $f(x+ne) = f(x) + nu$ and $f(x+(n+m)e) = f(x) + (n+m) u$ for some unit vector $u$. In particular, $f$ maps any isometric copy of ${\bf Z}$ to another isometric copy of ${\bf Z}$.

Next, if $e, f$ are orthogonal unit vectors in $X$, then $f(x), f(x+3e), f(x+4f)$ form a triangle with side lengths $3,4,5$ and in particular $f(x+3e)-f(x)$ is orthogonal to $f(x+4f)-f(x)$. From this and the previous claim we see that $f$ maps any isometric copy of ${\bf Z}^2$ to another isometric copy of ${\bf Z}^2$.

By the preceding discussion we have $$f(x,y)-f(0,0) =x (f(1,0)-f(0,0)) + y (f(0,1)-f(0,0)) \qquad (1)$$ whenever $x,y$ are integers, but also $f(n e)-f(0,0) = n (f(e)-f(0,0))$ for any unit vector $n$. In particular, for any Pythagorean triple $a^2+b^2=c^2$, we have (1) when $(x,y)$ is an integer multiple of $(\frac{a}{c},\frac{b}{c})$; conjugating by translation we then see that (1) also holds when $(x,y)$ is an integer multiple of $(\frac{a}{c},\frac{b}{c})$ plus an element of ${\bf Z}^2$. Iterating this we see that (1) holds whenever $(x,y)$ is an integer linear combination of Pythagorean fractions $(\frac{a}{c}, \frac{b}{c})$. In particular, (1) holds on a dense subset $D$ of ${\bf R}^2$.

This is already enough to get isometry in the continuous case. Now we remove the continuity hypothesis. Observe that if $\|x-y\| \leq 2$, then $x$ can be reached from $y$ by a pair of moves of unit length, and hence by the triangle inequality $\|f(x)-f(y) \| \leq 2$.

Now let us normalise $f,X,Y$ so that ${\bf R}^2 \subset Y$, $f(0,0)= (0,0)$, $f(1,0)=(1,0)$, and $f(0,1)=(0,1)$. Then by (1) we see that $f$ is the identity on $D$. For any $p \in {\bf R}^2$ and $\varepsilon>0$, we can find elements $q, r$ of $D$ within $\varepsilon$ of $p + (2,0)$ and $p-(2,0)$ that lie within $2$ of $p$, and thus by the preceding we see that $f(p)$ lies within $2+\varepsilon$ of $p+(2,0)$ and $p-(2,0)$. Sending $\varepsilon$ to zero we see that $f(p)=p$, and the claim follows.