EDIT 2: Made explicit the underlying non-asymptotic bounds in the calculation.
EDIT: Replaced the calculation for two dimensions by the case of arbitrary constant dimension. Added a table of the values.
I'd like to add an informal sketch of how David's very elegant result can be calculated. (To be clear, I suggest selecting his answer as the "correct" answer; this one is just intended to complement his.)
Assume the points are in general position, so that no two distinct pairs have the same distance. This happens with probability 1.
In the directed nearest-neighbor graph, each point has out-degree 1 (by definition). Also, for any directed path $p_1 \rightarrow p_2 \rightarrow p_3 \rightarrow \cdots \rightarrow p_k$ with no 2-cycles, we have $d(p_1, p_2) > d(p_2, p_3) > \cdots > d(p_{k-1}, p_k)$. That is, edge lengths decrease along the path. This is because, e.g., $p_3$ must be closer to $p_2$ than $p_1$ is (otherwise $p_3$ would not be $p_2$'s nearest neighbor), and so on.
As a consequence, in the undirected multigraph obtained by replacing each directed edge by its undirected equivalent, the only cycles are 2-cycles, where points $p_i$ and $p_j$ form a 2-cycle if and only if they are mutual nearest neighbors. Other edges are not on cycles.
It follows that the undirected nearest-neighbor graph (in which each such 2-cycle is replaced by a single edge) is acyclic, and has a number of edges equal to the number of vertices minus the number of pairs that are mutual nearest neighbors. Thus the number of components is equal to the number of pairs that are mutual nearest neighbors.
This holds in any metric space. Next, for intuition, consider the case of points in $R^1$, where calculations are relatively easy.
One Dimension
To ease the calculation, modify the distance metric to "wrap around" the boundary, that is, use
$$d_1(x, x') = \min\{|x-x'|,1-|x-x'|\}.$$ This changes the number of nearest-neighbor pairs by at most 1.
We need to estimate the expected number of pairs, among the $n$ points, that are mutual nearest neighbors. If we order the points as
$$p_1 < p_2 < \cdots < p_{n},$$
only pairs of the form $(p_i, p_{i+1})$ (or $(p_n, p_1)$) can be nearest neighbors. A given pair of this form are nearest-neighbors if and only if their distance is less than the distances of the neighboring pairs $(p_{i-1}, p_i)$ and $(p_{i+1}, p_{i+2})$ (to the left and right). That is, if, among the three consecutive pairs, the middle pair is closest. By symmetry (?), this happens with probability 1/3. Hence, by linearity of expectation, the number of the $n$ adjacent pairs that are nearest-neighbor pairs is $n/3$ (plus or minus 1, to correct for the wrap-around assumption). Hence the number of components is $n/3\pm 1$.
The symmetry argument above is suspect -- maybe there is some conditioning? It also doesn't extend to higher dimensions. Here's a more careful, detailed calculation that addresses these issues.
Let $p_1, p_2, \ldots, p_n$ be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs $n\choose 2$ times the probability that a given random pair $(p, q)$ are a nearest-neighbor pair. WLOG we can assume that $p$ and $q$ are the first two points drawn. Let $d_{pq}$ be their distance. They will be a nearest-neighbor pair if and only if none of the $n-2$ subsequent points are within distance $d_{pq}$ of $p$ or $q$.
The probability of this event (conditioned on $d_{pq}$) is $$\max(0, 1-3d_{pq})^{n-2},$$
because it happens if and only if none of the remaining $n-2$ points falls between $p$ and $q$ or within either of the two $d_{pq}$-wide boundaries on each side of $p$ and $q$.
To distance $d_{pq}$ is distributed uniformly over $[0, 1/2]$ (using our "wrap-around" assumption). Hence, the probability that $(p,q)$ is a nearest-neighbor pair is
$$\int_{0}^{1/3} (1-3 z)^{n-2} 2dz.$$
By a change of variables $x = 1-3z$ this is
$$\int_{0}^1 x^{n-2} 2\,dx/3 = \frac{2}{3(n-1)}.$$
By linearity of expectation, the expected number of nearest-neighbor pairs is $2{n\choose 2}/(3(n-1)) = n/3$ (plus or minus 1 to correct for the wrap-around technicality). Therefore the expected number of components is indeed $n/3\pm 1$.
As an aside, note that when $d_{pq}$ is large (larger than $\log(n)/n$, say), the contribution to the expectation above is negligible. So, we could under- or over-estimate the conditional probability for such $d_{pq}$ significantly; this would change the result by lower-order terms.
Any constant dimension
Fix constant dimension $k \in \{1,2,\ldots\}$.
To ease the calculations, modify the distance metric to wrap around the borders, that is, use
$d_k(p, q) = \sqrt{\sum_{i=1}^k d_1(p_i, q_i)^2}$
for $d_1$ as defined previously. This changes the answer by at most an additive $o(n)$ with high probability and in expectation.
Define $\beta_k, \mu_k\in \mathbb R$ such that $\beta_k r^k$ and $\mu_k r^k$ are the volumes of, respectively, a ball of radius $r$ and the union of two overlapping balls of radius $r$ whose centers are $r$ apart (so each center lies on the boundary of the other ball).
Let $p_1, p_2, \ldots, p_n$ be the points in sampled order. By linearity of expectation, the expected number of nearest-neighbor pairs is the number of pairs $n\choose 2$ times the probability that a given random pair $(p, q)$ are a nearest-neighbor pair. WLOG we can assume that $p$ and $q$ are the first two points drawn. Let $d_{pq}$ be their distance. They will be a nearest-neighbor pair if and only if none of the $n-2$ subsequent points are within distance $d_{pq}$ of $p$ or $q$.
We calculate the probability of this event. In the case that $d_{pq} \ge 1/4$, the probability of the event is at most the probability that no point falls within the ball of radius $1/4$ around $p$, which is at most $(1-\beta_k/4^k)^{n-2} \le \exp(-(n-2)\beta_k/4^k)$.
The case that $d_{pq} \le 1/4$ happens with probability $\beta_k/4^k$. Condition on any such $d_{pq}$. Then $p$ and $q$ will be nearest neighbors iff none of the $n-2$ subsequent points fall in the "forbidden" region consisting of the union of the two balls of radius $d_{pq}$ with centers at $p$ and $q$. The area of this region is $\mu_k d_{pq}^k$ by definition of $\mu_k$ (using here that $d_{pq}\le 1/4$ and the metric wraps around), so the probability of the event in question is $(1-\mu_k d_{pq}^k)^{n-2}$.
Conditioned on $d_{pq} \in [0,1/4]$, the probability density function of $d_{pq}$ is $f(r) = k 4^k r^{k-1}$ (note $\int_{0}^{1/4} k 4^k r^{k-1} = 1$). Hence, the overall (unconditioned) probability of the event is
$$\frac{\beta_k}{4^k} \int_{0}^{1/4} (1-\mu_k r^k)^{n-2} k 4^kr^{k-1} \, dr ~+~ \epsilon(n,k)$$
where $$0 \le \epsilon(n, k) \le \exp(-(n-2)\beta_k /4^k).$$
Using a change of variables $z^k=1-\mu_k r^k$ to calculate the integral, this is
$$\frac{k \beta_k}{\mu_k} \int_{\alpha}^1 z^{k(n-1)-1} \, dz ~+~ \epsilon'(n, k)
= \frac{\beta_k}{\mu_k}\frac{1 + \epsilon'(n, k)}{n-1}$$
for constant $\alpha=(1-\mu_k/4^k)^{1/k}<1$
and "error term" $\epsilon'(n, k)$ satisfying
$$-\exp(-(n-1)\mu_k/4^k) ~\le~ \epsilon'(n, k) ~\le~ \exp(-(n-2)\beta_k /4^k)(n-1)\mu_k/\beta_k$$
so $\epsilon'(n, k) \rightarrow 0$ as $n\rightarrow\infty$.
By linearity of expectation, the expected number of nearest-neighbor pairs is
$$\frac{n\choose 2}{n-1}\frac{\beta_k}{\mu_k}(1+ \epsilon'(n,k))
= \frac{\beta_k}{2\mu_k}(1 + \epsilon'(n,k)) n,$$
where $\epsilon'(n, k) \rightarrow 0$ as $n\rightarrow\infty$.
Correcting for the wrap-around assumption adds a $\pm o(n)$ term.
Hence, asymptotically, the expected number of mutual nearest-neighbor pairs is $n\beta_k/(2\mu_k) + o(n)$.
Next we give more explicit forms for $\beta_k/(2\mu_k)$.
According to this Wikipedia entry,
$$\beta_k = \frac{\pi^{k/2}}{\Gamma(k/2+1)} \sim \frac{1}{\sqrt{\pi k}}\Big(\frac{2\pi e}{k}\Big)^{k/2}$$
where $\Gamma$ is Euler's gamma function,
with $\Gamma(k/2+1) \sim \sqrt{\pi k}(k/(2e))^{k/2}$ (see here).
Following the definition of $\mu_k$, the volume of the union of the two balls is twice the volume of one ball with a "cap" removed (where the cap contains the points in the ball that are closer to the other ball).
Using this math.se answer (taking $d=r_1=r_2=r$, so $c_1=a=r/2$) to get the volume of the cap, this gives
$$\mu_k = \beta_k (2 - I_{3/4}((k+1)/2, 1/2)),$$
where $I$ is the "regularized incomplete beta function". Hence, the desired ratio is
$$\frac{\beta_k}{2\mu_k} = \Big(4-2I_{3/4}\Big(\frac{k+1}2, \frac{1}{2}\Big)\Big)^{-1}.$$
I've appended below the first 20 values, according to WolframAlpha.
$$
\begin{array}{cc}
\begin{array}{|rcl|}
k & \beta_k / (2\mu_k) & \approx \\ \hline
1 & \displaystyle\frac{1}{3} & 0.333333 \\
2 &\displaystyle\frac{3 \pi}{3 \sqrt{3}+8 \pi} & 0.310752 \\
3 &\displaystyle\frac{8}{27} & 0.296296 \\
4 &\displaystyle\frac{6 \pi}{9 \sqrt{3}+16 \pi} & 0.286233 \\
5 &\displaystyle\frac{128}{459} & 0.278867 \\
6 &\displaystyle\frac{15 \pi}{27 \sqrt{3}+40 \pi} & 0.273294 \\
7 &\displaystyle\frac{1024}{3807} & 0.268978 \\
8 &\displaystyle\frac{420 \pi}{837 \sqrt{3}+1120 \pi} & 0.265577 \\
9 &\displaystyle\frac{32768}{124659} & 0.262861 \\
10 &\displaystyle\frac{420 \pi}{891 \sqrt{3}+1120 \pi} & 0.26067 \\ \hline
\end{array} &
\begin{array}{|rcl|}
k & \beta_k / (2\mu_k) & \approx \\ \hline
11 &\displaystyle\frac{262144}{1012581} & 0.258887 \\
12 &\displaystyle\frac{330 \pi}{729 \sqrt{3}+880 \pi} & 0.257427 \\
13 &\displaystyle\frac{4194304}{16369695} & 0.256224 \\
14 &\displaystyle\frac{5460 \pi}{12393 \sqrt{3}+14560 \pi} & 0.255228 \\
15 &\displaystyle\frac{33554432}{131895783} & 0.254401 \\
16 &\displaystyle\frac{120120 \pi}{277749 \sqrt{3}+320320 \pi} & 0.253712 \\
17 &\displaystyle\frac{2147483648}{8483550147} & 0.253135 \\
18 &\displaystyle\frac{2042040 \pi}{4782969 \sqrt{3}+5445440 \pi} & 0.252652 \\
19 &\displaystyle\frac{17179869184}{68107648041} & 0.252246 \\
20 & \displaystyle\frac{38798760 \pi}{91703097 \sqrt{3}+103463360 \pi} & 0.251904 \\ \hline
\end{array}
\end{array}
$$
Notably, for $k=20$ (and larger, in fact), WolframAlpha reports numerical values approaching 0.25, in contrast to the experimental results reported by OP which are much lower. Where is this discrepancy coming from?
