It turns out the identities I needed for resolving
$\int_{-1}^1P_n(x)T_j(x)\mathrm{d}x$
into a closed form was well-hidden in Abramowitz and Stegun and Gradshteyn and Ryzhik.
As I had mentioned in the edit to my original question, the integral is 0 if $j<n$ by virtue of the orthogonality property of the Legendre polynomials.
I now considered the following integral:
$\int_{-1}^1P_n(x)T_{n+k}(x)\mathrm{d}x\quad k\geq0$
To dispose of an elementary case first, I noted that $P_n(x)T_{n+k}(x)$ is an odd function iff $k$ is odd and even iff $k$ is even; the integral is therefore 0 for odd $k$.
The even $k$ case I had solved by making use of two identities: this series representation for $T_{n}(x)$ (also in Abramowitz and Stegun as 22.3.6), and an integral I derived from a more general form in Gradshteyn and Ryzhik:
$\int_0^1x^{n+2\rho}P_n(x)\mathrm{d}x=\frac{\left(2\rho+1\right)_n}{2^{n+1}\left(\rho+\frac1{2}\right)_{n+1}}$
where $\left(a\right)_n$ is the Pochhammer symbol. (The identity actually listed in G&R was an integral for a Gegenbauer (ultraspherical) polynomial, of which the Legendre polynomial is a special case.)
I only needed to retain terms in the series greater than or equal to $n$, again due to orthogonality of the Legendre polynomial. Applying the integral formula to each term (with an additional factor of 2 because the integrand is even), and feeding the resulting sum to Mathematica netted the following closed form:
$\int_{-1}^1P_n(x)T_{n+2k}(x)\mathrm{d}x=-\frac1{4}\frac{\left(n+2k\right)\Gamma\left(n+k\right)\Gamma\left(k-\frac1{2}\right)}{\Gamma\left(k+1\right)\Gamma\left(n+k+\frac{3}{2}\right)}$
(The original result returned by Mathematica 5.2 had nasty cosecant factors, which I disposed of using the reflection formula for the gamma function).
This can then be applied to the original integral with the three polynomials by exploiting the product-sum identity for the Chebyshev polynomial.
Regarding which point of view is preferable, I don't see that there's much difference between them. The Chebyshev polynomials map $[-1,1]$ to $[-1,1]$, the spread polynomials map $[0,1]$ to $[0,1]$, and they are conjugate under a linear map between $[-1,1]$ and $[0,1]$, so all their properties translate easily between the two frameworks.
However, I'd vote for Chebyshev polynomials as being somewhat more fundamental, due to orthogonality. The spread polynomials aren't orthogonal with respect to any measure, because they are nonnegative everywhere. To get orthogonality, one must subtract $1/2$, after which they become orthogonal with respect to $dx/\sqrt{x(1-x)}$ on the interval $[0,1]$. By contrast, the Chebyshev polynomials are already orthogonal with respect to $dx/\sqrt{1-x^2}$ on $[-1,1]$, with no subtraction needed. This isn't a big deal, since it just amounts to subtracting $1/2$, but it's nice not to have to do the subtraction.
Overall, there's nothing sacred about using the domain $[-1,1]$ for Chebyshev polynomials. Of course it aligns beautifully with trigonometry, but Chebyshev polynomials are important in many other settings (such as approximation theory) in which $[-1,1]$ plays no special role, and they are simply rescaled to fit the interval of interest. From that perspective, $[0,1]$ is just as good a domain. On the other hand, I see no gain from making the range $[0,1]$ as well, and one has to undo it to recover orthogonality.
Comments added in edit:
As for the factorizations, this amounts to factoring $T_n(x)$ (for Chebyshev polynomials) or $T_n(x)+1$ (for spread polynomials - not quite, see comments below). Both are interesting, since both the roots and the extrema of the Chebyshev polynomials are important.
In fact, $T_{2n}(x)= T_2(T_n(x))$ and hence $T_{2n}(x)+1 = 2T_n(x)^2$, so factoring spread polynomials includes factoring Chebyshev polynomials as the even-index case. (In the odd-index case, $T_{2n+1}(x)+1 = (T_{n+1}(x)+T_n(x))^2/(x+1)$, but I'm not certain how to interpret this.)
So I'd say factoring spread polynomials is more general but slightly more obscure. Definitely both are interesting, though.
There are combinatorial interpretations of the Chebyshev polynomials involving weighted monomer-dimer configurations (although the conditions are a little odd: see http://www.math.hmc.edu/~benjamin/papers/CombTrig.pdf). The analogous idea doesn't work out as nicely for spread polynomials, but maybe some other approach is more appropriate. It's worth noting that the Chebyshev polynomials have somewhat simpler coefficients. For example, $T_8(x)=128x^8-256x^6+160x^4-32x^2+1$ while $S_8(x)=-16384x^8 + 65536x^7 - 106496x^6 + 90112x^5 - 42240x^4 + 10752x^3 - 1344x^2 + 64x$.
Best Answer
Both the Legendre and Chebyshev polynomials are particular cases of Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$. A general connection formula of the type $$P_n^{(\gamma,\delta)}(x)=\sum_{k=0}^nc_{n,k}^{\gamma,\delta;\alpha,\beta}P_k^{(\alpha,\beta)}(x)$$ can be found on page 256 of the book [Mourad E.H. Ismail, Classical and quantum orthogonal polynomials in one variable, Encyclopedia of Mathematics and its Applications 98, Cambridge University Press, Cambridge, 2005].