Short version
Consider the two following observations:
Given a state $\rho$, the problem of finding purifications of $\rho$ is equivalent to that of finding matrices $A$ such that $\rho=AA^\dagger$. The purifications of $\rho$ are then the vectorisations of these $A$ (i.e. a vector $\Psi$ is a purification of $\rho$ iff its matrix of coefficients $\Psi_{ij}$ satisfies $\Psi\Psi^\dagger=\rho$).
In this notation, the overlap between two purifications $A,B$ is written as $\operatorname{Tr}(A^\dagger B$).
We want to find the value of $\operatorname{Tr}|\sqrt\rho\sqrt\sigma|$.
As shown in this other answer, this is equivalent to finding the unitary $U$ maximising $\lvert\operatorname{Tr}(U\sqrt\rho\sqrt\sigma)\rvert$, and the unitary achieving this maximum turns out to equal $V^\dagger$ where $V$ is the unitary in the polar decomposition of $\sqrt\rho\sqrt\sigma$.
Observation (2) above tells us why finding the fidelity reduces to a maximisation problem, and how the polar decomposition (or more generally, the singular values) enters into it, while observation (1) tells us why we can understand the terms in this maximisation as overlaps of purifications of the states.
Consider then the $U$ such that $F(\rho,\sigma)=\operatorname{Tr}(U\sqrt\rho\sqrt\sigma)$, and define $A\equiv U\sqrt\rho$ and $B=\sqrt\sigma$. Then $F(\rho,\sigma)=\operatorname{Tr}(AB^\dagger)$, and thus $F(\rho,\sigma)$ equals the overlap between the purification of $\rho$ corresponding to $A$ and the purification of $\sigma$ corresponding to $B$.
Detailed version
The starting observation is that there are two "natural" ways to write down $\sqrt\rho\sqrt\sigma$ (or any other product of two normal operators): using their spectral decompositions, and using the singular value decomposition of their product:
$$\sqrt\rho\sqrt\sigma=\sum_{jk} \sqrt{\lambda_j\mu_k}\lvert\lambda_j\rangle\!\langle\lambda_j\rvert \mu_k\rangle\!\langle \mu_k\rvert
= \sum_m s_m \lvert s_m^L\rangle\!\langle s_m^R\rvert,$$
where $\rho=\sum_j\lambda_j\lvert\lambda_j\rangle\!\langle\lambda_j\rvert$ and $\sigma=\sum_k\mu_k\lvert\mu_k\rangle\!\langle\mu_k\rvert$ are the spectral decompositions of $\rho$ and $\sigma$, and I denoted with $s_m$ the singular values of $\sqrt\rho\sqrt\sigma$, and with $\lvert s_m^{L(R)}\rangle$ the left (right) singular vectors of $\sqrt\rho\sqrt\sigma$.
Note that using these definitions, we have
$\lvert\sqrt\rho\sqrt\sigma\rvert=\sum_m s_m \lvert s_m^R\rangle\!\langle s_m^R\rvert$ (if using the definition $\lvert A\rvert\equiv\sqrt{A^\dagger A}$, otherwise replace $R$ with $L$ if you want to define $\lvert A\rvert\equiv\sqrt{AA^\dagger}$),
and thus $\operatorname{tr}\lvert\sqrt\rho\sqrt\sigma\rvert=\sum_m s_m$.
Let us now denote with $\lvert\psi_\rho\rangle$ and $\lvert\psi_\sigma\rangle$ a pair of purifications of $\rho$ and $\sigma$. These can be written in general as
$$\lvert\psi_\rho\rangle=\sum_k \sqrt{\lambda_k}\lvert\lambda_k\rangle\otimes\lvert u_k\rangle, \\
\lvert\psi_\sigma\rangle=\sum_k \sqrt{\mu_k}\lvert\mu_k\rangle\otimes\lvert v_k\rangle,
$$
for arbitrary orthonormal bases $\{u_k\}_k, \{v_k\}_k$. We can then write their overlap as
$$\langle\psi_\rho\rvert\psi_\sigma\rangle = \sum_{jk}\sqrt{\lambda_j\mu_k}
\langle\lambda_j\rvert\mu_k\rangle\langle u_k\rvert v_k\rangle =
\sum_{jk}\langle\lambda_j\rvert\sqrt\rho\sqrt\sigma\lvert\mu_k\rangle \langle u_j\rvert v_k\rangle,
$$
where I exploited the fact that $\sqrt\rho\lvert\lambda_j\rangle=\sqrt{\lambda_j}\lvert\lambda_j\rangle$ and
$\sqrt\sigma\lvert\mu_k\rangle=\sqrt{\mu_k}\lvert\mu_k\rangle$.
Using the singular value decomposition of $\sqrt\rho\sqrt\sigma$ we thus get
\begin{align}\langle\psi_\rho\rvert\psi_\sigma\rangle &= \sum_{jkm}
s_m \langle\lambda_j\rvert s_m^L\rangle\!\langle s_m^R\rvert\mu_k\rangle \langle u_j\rvert v_k\rangle \\
&= \sum_m s_m \langle s_m^R\rvert
\Bigg(\underbrace{\sum_k \lvert \mu_k\rangle\!\langle \bar{v}_k\rvert }_{U_1}\Bigg)
\Bigg(\underbrace{\sum_j \lvert\bar{u}_j\rangle\!\langle \lambda_j\rvert}_{U_2}\Bigg)
\lvert s_m^L\rangle,
\end{align}
where I denoted with $\lvert \bar{u}_j\rangle$ the complex conjugate vector of $\lvert u_j\rangle$, so that $\langle u_j\rvert v_k\rangle=\langle \bar{v}_k\rvert \bar{u}_j\rangle$.
Uhlmann's theorem is almost straightforward from this.
The triangle inequality gives
$\lvert\langle\psi_\rho\rvert\psi_\sigma\rangle\rvert\le \sum_m s_m$ because matrix elements of unitary matrices are always less than $1$ in modulus, and the inequality is saturated when
$U_1 U_2=\sum_m \lvert s_m^R\rangle\!\langle s_m^L\rvert\equiv \mathcal U_{PD}^\dagger$.
In terms of the purification vectors $\lvert u_j\rangle,\lvert v_k\rangle$, this happens when
$$\langle u_j\rvert v_k\rangle = \langle \mu_k \rvert \Bigg(\sum_m \lvert s_m^R\rangle \!\langle s_m^L\rvert \Bigg) \lvert \lambda_j\rangle
= \langle \mu_k\rvert \mathcal U_{PD}^\dagger\rvert \lambda_j\rangle.$$
Note that here $\mathcal U_{PD}$ is the unitary matrix that you get out of the polar decomposition, that is, the $V$ in your post.
We can thus conclude that the purifications that saturate the inequality are those of the form
\begin{align}
\lvert u_j\rangle = V\lvert\bar{\lambda}_j\rangle, \qquad
\lvert v_k\rangle = V\mathcal U_{PD}^*\lvert\bar{\mu}_k\rangle,
\end{align}
or, equivalently,
\begin{align}
\lvert u_j\rangle = V\lvert j\rangle, \qquad
\lvert v_k\rangle = V\lvert j\rangle \langle \mu_k\rvert\mathcal U_{PD}^\dagger\lvert\lambda_j\rangle,
\end{align}
for any unitary $V$.
So, in conclusion, what does this tell us? That the vectors $\lvert u_j\rangle,\lvert v_k\rangle$ that make the purifications the most aligned are determined by the overlap between the eigenvectors of $\sigma$, and the eigenvectors of $\rho$ rotated through unitary that maps the right singular vectors of $\sqrt\rho\sqrt\sigma$ into the left ones.
Why is this the case? I have no idea, mostly because I don't know of any easy way to relate the SVD of a product of two operators with their eigenvectors.
