The answer turns out to be
$$ \frac{4\pi}{5} \approx 2.513274123. $$
Below is a simulation of $10^7$ trials using Mathematica 13, demonstrating a sample mean that is quite close to the exact value above:
![Simulation](https://cdn.statically.io/img/i.sstatic.net/TJkgb.png)
Proof. Let $A, B, C$ be independent, uniformly sampled points on the unit sphere $\mathbb{S}^2 $ centered at the origin.
Step 1. It is easy to check that the distance $L$ between the origin and the line $\overline{AB}$ has the PDF of the form
$$ f_L(l) = \begin{cases}
2l, & 0 < l < 1, \\
0, & \text{elsewhere}.
\end{cases} $$
Indeed, assuming WLOG that $A = (0, 0, 1)$ is the north pole,
$$ \mathbf{P}(L \leq l) = \mathbf{P}( \text{[$z$-coordinate of $B$]} \leq 2l^2 - 1) = l^2 $$
and hence the desired claim follows.
Step 2. Now we condition on the event $\{L = l\}$. In doing so, we may assume that the line $\overline{AB}$ lies in the $xy$-plane, passes through the point $(l, 0, 0)$, and is perpendicular to the $x$-axis. In other words, we may assume that
$$ A = (l, \sqrt{1-l^2}, 0) \qquad\text{and}\qquad B = (l, -\sqrt{1-l^2}, 0). $$
Now, given the value of $l$, we introduce the following parametrization of $\mathbb{S}^2$:
$$ \Psi_l(r, \theta) = r \left(\sin\alpha, 0, \cos\alpha\right) + \sqrt{1-r^2} \bigl[ (\cos\alpha, 0, -\sin\alpha) \cos\theta + (0, 1, 0) \sin\theta \bigr] $$
Here, $-l \leq r \leq l$ and $-\pi \leq \theta \leq \pi$, and $\alpha = \arcsin(r/l)$. This parametrization is designed so that the following property holds:
Property. For each fixed $r$, the curve $\theta \mapsto \Psi_l(r, \theta)$ is a circle that passes through the points $A$ and $B$, and the plane containing this circle is at a distance of $|r|$ from the origin.
Below is an animated example of this circle where $l = 2/3$ and $r$ varies from $-l$ to $l$:
![Animated example](https://cdn.statically.io/img/i.sstatic.net/Z5czV.gif)
Then a bit of tedious computation shows that the area element is computed as
$$ \mathrm{d}A
= \left\| \frac{\partial \Psi_l}{\partial r} \times \frac{\partial \Psi_l}{\partial \theta} \right\| \, \mathrm{d}r\mathrm{d}\theta
= \left|1 - \sqrt{\frac{1-r^2}{l^2-r^2}}\cos\theta\right| \, \mathrm{d}r\mathrm{d}\theta. $$
So, if $R$ denotes the distance between the origin and the plane spanned by the points $A, B, C$, then
\begin{align*}
\mathbf{P}(R \leq \rho \mid L = l)
&= \frac{1}{4\pi} \int_{\operatorname{Dom}(\Phi_l)} \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}A \\
&= \frac{1}{4\pi} \int_{-l}^{l} \int_{-\pi}^{\pi} \left|1 - \sqrt{\frac{1-r^2}{l^2-r^2}}\cos\theta\right| \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}\theta\mathrm{d}r.
\end{align*}
Computing the inner integral by invoking the identity
$$ \int_{-\pi}^{\pi} |1 - k \cos\theta| \, \mathrm{d}\theta = 4\left( \frac{\pi}{2} - \arccos\left(\frac{1}{k}\right) + \sqrt{k^2 - 1} \right), $$
which holds for $k \geq 1$, the conditional CDF is recast as
\begin{align*}
\mathbf{P}(R \leq \rho \mid L = l)
&= \int_{-l}^{l} \frac{1}{\pi} \left( \frac{\pi}{2} - \arccos\sqrt{\frac{l^2-r^2}{1-r^2}} + \sqrt{\frac{1-l^2}{l^2-r^2}} \right) \mathbf{1}_{\{|r| \leq \rho\}} \, \mathrm{d}r.
\end{align*}
From this, the conditional PDF of $R$ given $L$ is computed as
$$ f_{R \mid L} (r \mid l) = \frac{2}{\pi} \left( \frac{\pi}{2} - \arccos\sqrt{\frac{l^2-r^2}{1-r^2}} + \sqrt{\frac{1-l^2}{l^2-r^2}} \right) \mathbf{1}_{\{r \leq l\}}. $$
Then by a tedious calculation (combined with integration by parts), we can check that the distribution of $R$ is fairly elementary, with the PDF given by
$$ f_{R}(r) = \frac{3}{2}( 1 - r^2 ). $$
Therefore the desired answer is
$$ \mathbf{E}[\pi(1-R^2)] = \frac{3\pi}{2} \int_{0}^{1} (1 - r^2)^2 \, \mathrm{d}r = \boxed{\frac{4\pi}{5}} \approx 2.513274123. $$