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
OP is essentially pondering if commutativity is a transitive relation, ie. if three normal operators$^{1}$ $A$, $B$, and $C$ satisfies
$$ [A,B]~=~0\quad \wedge\quad [B,C]~=~0 \quad\stackrel{?}{\Rightarrow} \quad [A,C]~=~0 .\tag{T}$$
The answer is No, but OP argues via the existence of a common basis of eigenvectors for two commuting normal operators that eq. (T) should hold.
To most clearly expose the flaw in OP's argument, pick $B$ to be proportional to the identity. Then $B$ commutes with everything. Clearly we can then find two non-commuting normal operators $A$ and $C$, so that eq. (T) is violated! And a basis of eigenvectors for $A$ cannot be a basis of eigenvectors for $C$, and vice versa.
--
$^{1}$We will ignore subtleties with unbounded operators, domains, selfadjoint extensions, etc., in this answer.