The integrable spin-boson model is not an inverse-square Ising model. For that equivalence to hold, one has to include terms that do not commute with X⊗1.
Let's re-examine the Feynman-Kac representation for the integral kernel of the partition function: when [H,X]=0, Kronecker deltas appear which massively simplify the partition function:
Z[τc] ∼τc→0∑{X(τ)},{xu(τ)}⋯⟨X(τi−1),xu(τi−1)|e−τc⋅X(τi)∑uλuxu(τi)e−τcHB|X(τi),xu(τi)⟩⋯
=∑{X(τ)},{xu(τ)}⋯δ(X(τi−1),X(τi))⋅⟨xu(τi−1)|e−τc⋅X(τi)∑uλuxu(τi)e−τcHB|xu(τi)⟩⋯
To justify the intermediate step, since HB does not affect the system degrees of freedom, X(τi) is just an eigenvalue, and so the bra ⟨X(τi−1)| meets the ket |X(τi)⟩ directly, inducing a prefactor. These Kronecker deltas in the integral kernel are precisely the reason why the thermal Green's function ⟨X(τ)X(τ′)⟩β for the integrable spin-boson model is identically one (as you correctly deduced by a direct calculation sans path integral techniques).
Of course, if we turned on a nontrivial system Hamiltonian HS (e.g. add a term proportional to Z⊗1), then there would be a possibility to connect different eigenspaces of X at different imaginary times, and these Kronecker-deltas would disappear (this is the starting point of the paper which you cited). But then the integrability of the model would be lost. Either way, the model is either solvable via a polaron transformation, or is equivalent to an inverse-square Ising model; both cannot be true simultaneously.