TL;DR: The key observation is that Schmidt basis on a subsystem consists of eigenvectors of the reduced state of that subsystem. Consequently, if the reduced state is a product state then its Schmidt basis can be chosen to consist of pure product states.
We begin by writing down the Schmidt decomposition for $|\psi\rangle_{ABC}$ with respect to the partitioning of $ABC$ into subsystems $B$ and $AC$
$$
|\psi\rangle_{ABC}=\sum_i s_i|\beta_i\rangle_B|\alpha_i\rangle_{AC}\tag1
$$
which allows us to write an eigendecomposition for $\rho_{AC}=\rho_A\otimes\rho_C$ in terms of the Schmidt basis
$$
\rho_A\otimes\rho_C=\sum_i s_i^2|\alpha_i\rangle_{AC}\langle \alpha_i|_{AC}.\tag2
$$
However, if
$$
\rho_A=\sum_k\lambda_k|\psi_k\rangle_A\langle\psi_k|_A\quad
\rho_C=\sum_l p_l|\varphi_l\rangle_C\langle\varphi_l|_C\tag3
$$
are eigendecompositions of $\rho_A$ and $\rho_C$, then
$$
\rho_A\otimes\rho_C=\sum_{kl}\lambda_k p_l|\psi_k\rangle_A|\varphi_l\rangle_C\langle\psi_k|_A\langle\varphi_l|_C\tag4
$$
is also an eigendecomposition of $\rho_A\otimes\rho_C$. Comparing $(2)$ and $(4)$, we see that both $s_i^2$ and $\lambda_k p_l$ are eigenvalues of $\rho_A\otimes\rho_C$. Therefore, the ranges of the index $i$ and the pair $k,l$ coincide, there is a bijection $f$ such that $f(k, l)=i$ and
$$
s_{f(k,l)}^2=\lambda_k p_l.\tag5
$$
Moreover, if the eigenvalues $s_i^2$ are all distinct, then
$$
|\alpha_{f(k,l)}\rangle_{AC}\equiv|\psi_k\rangle_A|\varphi_l\rangle_C\tag6
$$
where $\equiv$ denotes equality up to global phase. If there are repeated eigenvalues then $(6)$ does not necessarily hold. However, in this case, the unitary freedom$^1$ in the subspaces corresponding to equal Schmidt coefficients in $(1)$ allows us to choose $|\alpha_i\rangle_{AC}$ as in $(6)$.
Now, we can upgrade equality up to global phase in $(6)$ to proper equality by absorbing the phase factor into $|\beta_i\rangle_B$. Denoting the new basis, with phase factors absorbed, as $|\phi_{kl}\rangle_B$, we have
$$
\begin{align}
|\beta_{f(k,l)}\rangle_B|\alpha_{f(k,l)}\rangle_{AC}&=|\phi_{kl}\rangle_B|\psi_k\rangle_A|\varphi_l\rangle_C\\
&=|\psi_k\rangle_A|\phi_{kl}\rangle_B|\varphi_l\rangle_C.
\end{align}\tag7
$$
Finally, substituting $(5)$ and $(7)$ into $(1)$, we obtain
$$
\begin{align}
|\psi\rangle_{ABC}&=\sum_i s_i|\beta_i\rangle_B|\alpha_i\rangle_{AC}\\
&=\sum_{kl} s_{f(k,l)}|\beta_{f(k,l)}\rangle_B|\alpha_{f(k,l)}\rangle_{AC}\\
&=\sum_{kl}\sqrt{\lambda_k p_l}|\psi_k\rangle_A|\phi_{kl}\rangle_B|\varphi_l\rangle_C
\end{align}\tag8
$$
as desired.
$^1$ This is analogous to unitary freedom in singular value decomposition from which it derives.
