When you refer to SPTO transitions as being 'nonlocal symmetry-breaking transitions', I guess you are referring to the idea of 'hidden symmetry-breaking'? This indeed the old way of looking at SPTOs (i.e., 80s and 90s), before the more modern understanding in terms of topological invariants was developed. For clarity, let me briefly review this with an illustrative example before then addressing your questions:
Minimal SPT example in 1D
The minimal SPT phase in one dimension is the Haldane phase, e.g., protected by $\mathbb Z_2 \times \mathbb Z_2$ symmetry. A minimal exactly-solvable spin chain in this phase is the cluster model:
$$ H_\textrm{cluster} = -\sum_n X_{n-1} Z_n X_{n+1},$$
where I denote the spin-$1/2$ Pauli matrices as $X,Y,Z$. The protecting symmetries are $P_1 = \prod_n Z_{2n-1} $ and $P_2 = \prod_n Z_{2n}$. It turns out that the ground state has long-range order in a string order parameter:
$$ \lim_{|n-m|\to \infty} \langle X_{2m} Z_{2m+1} Z_{2m+3} \cdots Z_{2n-3} Z_{2n-1} X_{2n} \rangle \neq 0. $$
Note that the endpoint operator of the string (e.g., $X_{2n}$ on the right) is charged under the $P_2$ symmetry. Hence, it looks a bit like a symmetry-breaking order parameter, except for the nonlocal ($P_1$) string connecting the two endpoints.
''Hidden'' symmetry-breaking
Now, there is in fact a nonlocal mapping that makes this into a symmetry-breaking phase: define
$$ \tilde X_n = \cdots Z_{n-5} Z_{n-3} Z_{n-1} X_n \qquad \textrm{and} \qquad \tilde Z_n = Z_n.$$
Ignoring boundary issues, this gives a well-defined Pauli algebra. In these new variables, we have
$$ H_\textrm{cluster} = - \sum_n \tilde X_{n-1} \tilde X_{n+1}. $$
It is just an Ising chain on the even and odd sites! Moreover, the above string order parameter now just becomes
$$ \lim_{|n-m|\to \infty} \langle \tilde X_{2m} \tilde X_{2n} \rangle \neq 0, $$
i.e., it is just the conventional symmetry-breaking order parameter. If you wish, you can see this as a derivation of the earlier claim that the original cluster model has long-range order in that string order parameter (since in these effective Ising variables it is obvious that we have the local order parameter).
Pros and cons
A positive feature of this mapping is that it makes readily apparent that, e.g., the cluster model $H_\textrm{cluster}$ cannot be adiabatically connected to the trivial Hamiltonian $H_\textrm{triv} = - \sum_n Z_n$ (at least if we preserve $\mathbb Z_2 \times \mathbb Z_2$ symmetry). Indeed, under the above nonlocal mapping, the latter maps to $H_\textrm{triv} = - \sum_n \tilde Z_n$, which is symmetry-preserving and must thus be separated from the $-\sum_n \tilde X_{n-1} \tilde X_{n+1}$ phase by a thermodynamic phase transition.
The downside is that the nonlocal mapping hides the physics of the original model. After all, the cluster model has no symmetry-breaking, in particular, its ground state is unique for periodic boundary conditions. (This is possible due to the nonlocal mapping we used, which inevitably screws up when attempting to deal with periodic boundary conditions.) A more modern understanding of SPT phases is built on the concept of symmetry fractionalization (which I don't really want to go into here, if not just for space constraints, but for further reading I can point you to the intro of my paper here: p2 and p3 of arxiv:1707.05787; disclaimer: I did not invent the concept of symmetry fractionalization, and all the appropriate references can be found there).
Answering your questions
Is the symmetry breaking between SPTO phases always signalled by a non-local order parameter?
Roughly speaking: yes. More precisely, for 1D SPT phases this general construction was done in 2012 by Pollmann and Turner ( arxiv:1204.0704 ). This works for the vast majority of SPT phases, although there are funky exceptions discussed in the paper (which have thus far not appeared in any concrete model in the literature). Relatedly, Else, Barlett and Doherty ( arxiv/1304.0783 ) generalized the above nonlocal mapping from SPT to symmetry-breaking phases. Moreover, in higher dimensions the notion of symmetry fractionalization also allows to prove that one has long-range order for a symmetry membrane operator dressed with an appropriate unitary on its boundary (i.e., the higher-dimensional analogue of what we saw above for the cluster model, where the string had some endpoint operators). However, in that case I am not aware of a dual 'hidden symmetry breaking' picture (which is yet another downside of this way of understanding SPT phases).
Conversely, is symmetry breaking with a non-local order parameter always SPTO (there exists a gapless path...)
Sticking to the simpler one-dimensional case for convenience: basically, yes. More precisely, if you have an on-site symmetry $\prod_n U_n$ which has long-range order with some endpoint operator $\mathcal O_n$, and this endpoint operator is non-trivially charged under some (other) symmetry of your model, then you have a non-trivial SPT phase on your hands (the one caveat is that you have to make sure that this charge is discrete and cannot be smoothly connected to zero, e.g., $-1$ is not a discrete charge for $U(1)$, but it is a discrete charge for $\mathbb Z_2$).
In general, what is the relation between non-local order parameters
and topological phases?
I hope the above comments already address this point, at least for the case of SPT phases. Even intrinsic topological order can be related to the long-range order of certain string order parameters, known as the Bricmont-Frohlich-Fredenhagen-Marcu order parameter, which is conceptually completely unrelated to the above SPT case (for a condensed matter perspective, see this paper by Gregor, Huse, Moessner and Sondhi: arxiv:1011.4187 ).
