I also know that L and S commute, but I am unsure why. I've heard that it is simply because they act on difference variables, but I don't understand exactly what this means. Is there a way to show this explicitly?
Suppose we have two Hilbert spaces $H_1$ and $H_2$, an operator $A_1$ acting on $H_1$, and an operator $A_2$ acting on $H_2$. Let $H = H_1 \otimes H_2$. Then we can define $A_1$ and $A_2$ on $H$ by defining
\begin{align}
A_1(|a_1\rangle \otimes |a_2\rangle) &= A_1|a_1\rangle \otimes |a_2\rangle \\
A_2(|a_1\rangle \otimes |a_2\rangle) &= |a_1\rangle \otimes A_2|a_2\rangle
\end{align}
where $|a_1\rangle \in H_1, |a_2\rangle \in H_2$; and extending linearly to all of $H$. Then
\begin{align}
A_1 \circ A_2(|a_1\rangle \otimes |a_2\rangle) &= A_1(|a_1\rangle \otimes A_2|a_2\rangle) \\
&= A_1|a_1\rangle \otimes A_2|a_2\rangle \\
&= A_2(A_1|a_1\rangle \otimes |a_2\rangle) \\
&= A_2 \circ A_1(|a_1\rangle \otimes |a_2\rangle)
\end{align}
so the commutator vanishes on all pure tensors, and hence on all of $H$.
This is precisely the situation we have with the operators $L$ and $S$. In general, the wave function of a particle lives in a tensor product space. The spatial part of the wave function lives in one space, that of square-integrable functions on $\mathbb{R}^3$. The spin part, on the other hand, lives in a spinor space, i.e., some representation of $SU(2)$. $L$ acts on the spatial part, whereas $S$ acts on the spin part.
What are the remaining commutation relations between $J$, $J^2$, $L$, $L^2$, $S$, and $S^2$?
You should be able to work these out on your own, using the commutation and anti-commutation relations you already know, and properties of commutators and anti-commutators. For example,
$$[J_i, L_j] = [L_i + S_i, L_j] = [L_i, L_j] + [S_i, L_j] = i\hbar\epsilon_{ijk} L_k$$
Likewise:
\begin{align}
[J_i^2, L_j] &= J_i[J_i, L_j] + [J_i, L_j]J_i \\
&= J_i(i\hbar \epsilon_{ijk}L_k) + (i\hbar\epsilon_{ijk}L_k)J_i \\
&= i\hbar\epsilon_{ijk}\{J_i, L_k\} \\
&= i\hbar\epsilon_{ijk}\{L_i + S_i, L_k\} \\
&= i\hbar\epsilon_{ijk}(\{L_i, L_k\} + \{S_i, L_k\}) \\
&= i\hbar\epsilon_{ijk}(2\delta_{ik} I + 2S_i L_k) \\
&= 2i\hbar\epsilon_{ijk} S_i L_k
\end{align}
where we have used the fact that $L_i$ and $S_j$ commute, the linearity of $[,]$ and $\{,\}$, and the identity
$$[AB, C] = A[B, C] + [A, C]B$$
Best Answer
Since both the $L_i$ and the $A_i$ are constants of motion, their Poisson bracket with the Hamiltonian vanishes and so it is "as good as a constant" for purposes of the algebraic structure - like a constant, it commutes with everything. Essentially you get one differently scaled $\mathfrak{so}(4)$ algebra for each subspace defined by constant energy $H(q,p) = E$.