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$$
Vectors transform linearly,
\begin{equation}
x _i \rightarrow A _{ ij } x _j
\end{equation}
through some transformation matrix $A$.
Now consider the transformation of $ M _i $:
\begin{align}
e ^{ i \theta _j J _j } M _i e ^{ - i J _j \theta _j } & = \left( 1 + i \theta _j J _j - ... \right) M _i \left( 1 - i \theta _j J _j - ... \right) \\
& = M _i + i \theta _j \left[ J _j , M _i \right] + ...
\end{align}
Now if $ \left[ J _j , M _i \right] $ is only proportional to $ M _j $ as above then infinitesimally,
\begin{align}
e ^{ i \theta _j J _j } M _i e ^{ - i J _j \theta _j } & = M _i + i \theta _j \epsilon _{ jik} M _k \\
&= (\delta_{ik}+i\theta_j\epsilon_{jik})M_k\\
&= A _{i,j}M _j
\end{align}
for some transformation matrix $A$ as required.
Best Answer
Under a linear transformation $T$ of the vector space, operators $O$ on it transform as $O\mapsto TOT^\dagger$. Since by definition the $J_i$ are the infinitesimal generators of rotation as $R(\phi) = \mathrm{e}^{\mathrm{i}J_i\phi}$, the finite rotation $\mathrm{e}^{\mathrm{i}J_i \phi}O\mathrm{e}^{-\mathrm{i}J_i\phi}$ implies that the infinitesimal change of any observable under rotation is $[J_i,O]$. (This follows by Taylor expanding the exponentials and keeping only the term to first order in $\phi$.)
So the commutation relations of the form $[J_i,O]$ tell you how $O$ changes under rotation.
The exact same reasoning goes through for the infinitesimal generators $K_i$ of the Lorentz boosts, so $[K_i,H]$ is the infinitesimal change of energy under a boost.